Ice thickness from LiDAR

Estimating ice thickness from LiDAR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes on ice thickness estimation from airborne LiDAR\n",
    "Adam Steer, January 2016.\n",
    "adam.d.steer@gmail.com"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this project, we have a bunch of LIDAR points georeferenced relative to an ITRF08 ellipsoid. What we need is the height of each point relative to sea level.\n",
    "\n",
    "But where exactly is the sea surface? We can estimate it from a gravity model, plus a tide model, plus some knowledge of dynamic sea surface topography (wind pressure etc. ).\n",
    "\n",
    "Or we can use a sensor-based approach, since we know some points are returns from water, or extremely thin new ice.\n",
    "\n",
    "This document describes the second approach.\n",
    "\n",
    "Using a small subset of LiDAR as an example, this notebook steps through the process of estimating sea ice thickness."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, import some libraries we will need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#setup\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from datetime import datetime\n",
    "\n",
    "#Used for some functions involving finding points within\n",
    "# N metres of any other point\n",
    "from sklearn import neighbors as nb\n",
    "\n",
    "#Testing a Kriging library from sklearn. Memory hungry!\n",
    "#from sklearn.gaussian_process import GaussianProcess\n",
    "\n",
    "from scipy.interpolate import SmoothBivariateSpline\n",
    "#Scipy's least squares bivariate spline was tested\n",
    "# but less functional than it's higher level SmoothBivariateSpline\n",
    "#from scipy.interpolate import LSQBivariateSpline\n",
    "\n",
    "#used to find the N closest points to neighbourhood median, \n",
    "# for filtering 'water' keypoints\n",
    "from heapq import nsmallest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function declarations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this function is the snow depth model. Parameters come from Steer et al (2016):\n",
    "#Estimating snow depth from altimery for East Antarctic pack ice\n",
    "\n",
    "def compute_zs(tf, s_i, tf_uncert):\n",
    "    \"\"\"\n",
    "    Take in the total freeboard (tf, float array), slope and intercept for an empirical model\n",
    "    of snow depth from elevation (s_i, tuple) and the total freeboard uncertainty\n",
    "    (tf_uncert, float array)\n",
    "\n",
    "    Return a snow depth, with an associated uncertainty.\n",
    "    zs, uncert = compute_zs(tf, ([slope, intercept]), tf_uncert)\n",
    "    \"\"\"\n",
    "    \n",
    "    zs = (params[0] * tf) + s_i[1]\n",
    "\n",
    "    zs_uncert = s_i[0] * tf_uncert\n",
    "\n",
    "    return zs, zs_uncert\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# next define the model for estimating ice thickness from total freeboard, snow depth\n",
    "# and some density parameters:\n",
    "\n",
    "def compute_zi(tf, zs, d_ice, d_water, d_snow, sd_tf, sd_zs, sd_dsnow, \\\n",
    "               sd_dice, sd_dwater):\n",
    "    \"\"\"\"\n",
    "    sea ice thickness from elevation, and propagation of uncertainties\n",
    "    after Kwok, 2010; kwok and cunningham, 2008\n",
    "    equations:\n",
    "    4: sea ice thickness from elevation\n",
    "    6: taylor series expansion of variance/covariance propagation using\n",
    "    partial derivatives\n",
    "    8, 9 and 10: partial derivatives used in the taylor series\n",
    "    \"\"\"\n",
    "    zi = (d_water / (d_water-d_ice)) * tf - ((d_water-d_snow) / \\\n",
    "         (d_water - d_ice)) * zs\n",
    "\n",
    "    zi_uncert = sd_tf**2 * (d_water / (d_water - d_ice))**2 + \\\n",
    "           sd_zs**2 * ((d_snow - d_water) / (d_water - d_ice))**2 + \\\n",
    "           sd_dsnow**2 * (zs / (d_water - d_ice))**2 + \\\n",
    "           sd_dice**2 * (tf /  (d_water - d_ice))**2 + \\\n",
    "           sd_dwater**2 * (((-d_ice * tf) + ((d_ice-d_snow) * zs)) / \\\n",
    "           (d_water - d_ice)**2)**2\n",
    "\n",
    "    return zi, zi_uncert\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# A function to build a KD-tree\n",
    "def build_tree(pointcloud, leafsize):\n",
    "    return nb.KDTree(pointcloud, leaf_size=leafsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#part one of fiding water - intensity percentile based low pass filter\n",
    "def find_low_intens(intens,pntile):\n",
    "    \"\"\"\n",
    "    Take in some return intensity data (array) and a percentile (float)\n",
    "    Return the indices of chosen percentile of intensity values\n",
    "    Some of these will be returns on ice, which need to be exlcuded.\n",
    "    A number-of-neighbours test is used to do this.\n",
    "    \"\"\"\n",
    "    li_idx =  np.where(intens <= np.percentile(intens, pntile))\n",
    "    return list(li_idx[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# part two of finding sea level keypoints\n",
    "def find_lowpoints(xyz, leafsize, nhood):\n",
    "    \"\"\"\n",
    "    Pass in an array of elevation vlues and an integer  number of points to use\n",
    "    as a neighbourhood.\n",
    "    Get back the indices of 10 elevation values closest to the neighbourhooed median\n",
    "    for neighbourhoods with more than nhood*4 points\n",
    "    It's not really finding low points, so the function needs renaming.\n",
    "    \"\"\"    \n",
    "    \n",
    "    low_tree = build_tree(xyz, leafsize)\n",
    "    nhoods = low_tree.query_radius(xyz, r=nhood)\n",
    "    point_filt = []   \n",
    "  \n",
    "    for nh in nhoods:\n",
    "        #print(nhood)\n",
    "        #print(xyz[nhood,2])\n",
    "        # we want indices in the input cloud\n",
    "        # of the 10 values closest to the nhood median\n",
    "        if len(nh) > 3*nhood:\n",
    "            npoints = sorted(enumerate(xyz[nh[:],2]),\\\n",
    "                      key=lambda x:abs(np.median(xyz[nh[:],2])-x[1]))[:10]\n",
    "            #print(npoints)\n",
    "            for point in npoints:\n",
    "                #print(point[0])\n",
    "                point_filt.append(nh[point[0]])\n",
    "                #print(xyz[point[0],:])\n",
    "\n",
    "            #point_filt.append(nh[find_nearest(xyz[nh[:],2], np.median(xyz[nh[:],2]))])\n",
    "        #print(new_z[i])\n",
    "    #return a list of indices\n",
    "    return point_filt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# gluing it together, using the two functions above\n",
    "def find_water_keys(xyzi, e_p, i_p):\n",
    "    \"\"\"\n",
    "    Pass in:\n",
    "    - a set of X,Y,Z and Intensity points\n",
    "    - the number of points to use as a neighbourhood for choosing low elevation\n",
    "    - a percentile level for intensity thresholding\n",
    "    Get back:\n",
    "    - a set of XYZI points corresponding to 'local sea level'\n",
    "    \"\"\"\n",
    "    #find lowest intensity values\n",
    "    low_intens_inds = find_low_intens(xyzi[:, 3], i_p)\n",
    "    li_xyzi = xyzi[low_intens_inds, :]\n",
    "\n",
    "    #find low elevations\n",
    "    low_elev_inds = find_lowpoints(li_xyzi[:, 2], e_p)\n",
    "    low_xyzi = li_xyzi[low_elev_inds, :]\n",
    "    #low_xyzi = np.squeeze(low_xyzi)\n",
    "\n",
    "    return low_xyzi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#now we have water keys, this function fits a 2D spline, the adjusts the points by evaluating the spline\n",
    "# at keypoints and subtracting the result from LiDAR points. The derived splines are also evaluated at *every*\n",
    "# lidar point, and the result subtracted - so there are some careful considerations about extrapolation from\n",
    "# polynomial splines and the shape of data being fitted. In this function, X/Y coordinates are normalised such that\n",
    "# [0 < X < 1] and [0 < Y < 1] before spline fitting. This reduces (but not destroys) the likelihood that\n",
    "# extrapolation of values outside the keypoint range results in wildly spurious values. As a precondition,\n",
    "# input data are clipped such that the along-track dimension is always within the range of keypoints\n",
    "# (ie swath chunks start and end over water).\n",
    "\n",
    "def f_points(in_, params, dims, smth):\n",
    "    '''\n",
    "    find open water points using find_water_keys\n",
    "    fit a surface to central tendency of keypoints\n",
    "    adjust the rest of a LiDAR swath by water keypoint values\n",
    "    In this method, coordinates are normalised so that \n",
    "    interpolation occurs in an almost-gridded environment which data pretty\n",
    "    much fill.\n",
    "    \n",
    "    Best applied with a lot of smoothing\n",
    "    '''\n",
    "    #normalise input data to a unit block\n",
    "    #http://stackoverflow.com/questions/3864899/resampling-irregularly-spaced-data-to-a-regular-grid-in-python\n",
    "    # first get min/max of points\n",
    "    xmin = np.min(in_[:,0])\n",
    "    ymin = np.min(in_[:,1])    \n",
    "    \n",
    "    #translate to starting at [0,0]\n",
    "    t_x = in_[:,0]-xmin\n",
    "    t_y = in_[:,1]-ymin\n",
    "\n",
    "    #normalise coordinates in each direction    \n",
    "    xmax = np.max(t_x)\n",
    "    ymax = np.max(t_y)\n",
    "    norm_x = t_x/xmax\n",
    "    norm_y = t_y/ymax\n",
    "    \n",
    "    #set up a translated world-space array to send to water point finding\n",
    "    to_waterkeys = np.column_stack([t_x, t_y, in_[:,2], in_[:,3]])\n",
    "    #print(to_waterkeys)\n",
    "    \n",
    "    #find water points   \n",
    "    xyz_fp = np.squeeze(find_water_keys(to_waterkeys, params[0], params[1]))\n",
    "    \n",
    "    #translate the results to the same normalised coordinate system\n",
    "    norm_fit_x = xyz_fp[:,0]/xmax\n",
    "    norm_fit_y = xyz_fp[:,1]/ymax\n",
    "    \n",
    "    #check extents\n",
    "    print('min. norm fit Y: {}, max. norm fit Y: {}'.format(min(norm_fit_y), max(norm_fit_y)))    \n",
    "    \n",
    "    #another 2D spline, pretty much the same as SmoothBivariateSpline\n",
    "    #xcoord = np.linspace(0, 1, 0.1)\n",
    "    #ycoord = np.linspace(0, 1, 0.1)\n",
    "    #this_f = LSQBivariateSpline(norm_fit_x, norm_fit_y,\n",
    "    #                            xyz_fp[:, 2], xcoord, ycoord)\n",
    "    \n",
    "    #this one is the winner right now.\n",
    "    #fit a spline using normalised XY coordinates\n",
    "    this_f = SmoothBivariateSpline(norm_fit_x, norm_fit_y,\n",
    "                                   xyz_fp[:, 2], kx=dims[0], ky=dims[1],\n",
    "                                   bbox=[0,1,0,1], s = smth*len(norm_fit_x))\n",
    "                                   \n",
    "    #evaluate the function at real-space coordinates\n",
    "    #fpoints_mod = this_f.ev(xyz_fp[:, 0], xyz_fp[:, 1])\n",
    "    #or normalised? Normalised, of course!\n",
    "    fpoints_mod = this_f.ev(norm_fit_x, norm_fit_y)\n",
    "    \n",
    "    #subtract evaluated spline heights from the LiDAR points\n",
    "    adjusted_points =  xyz_fp[:, 2] - fpoints_mod\n",
    "    \n",
    "    #one more filter - kill extreme values (+- 3 * sd)\n",
    "    e_f = np.where((adjusted_points >= 0-3*np.std(adjusted_points)) & (adjusted_points<= 0+3*np.std(adjusted_points)))\n",
    "\n",
    "    print('mean fitpoint adjustment (Z): {}'.format(np.mean(fpoints_mod[e_f[0]])))\n",
    "    \n",
    "    #translate fit points back!\n",
    "    fitpoints_mod = np.column_stack([xyz_fp[e_f[0], 0]+xmin, xyz_fp[e_f[0], 1]+ymin,\n",
    "                                    adjusted_points[e_f[0]]])\n",
    "    #evaluate the surface at \n",
    "    #z_mod = this_f.ev(in_[:, 0], in_[:, 1])\n",
    "    \n",
    "    #normalised coordinates\n",
    "    z_mod = this_f.ev(norm_x, norm_y)\n",
    "    \n",
    "    coeffs = this_f.get_coeffs()\n",
    "    resid = this_f.get_residual()    \n",
    "    \n",
    "    #remove predicted values from elevations\n",
    "    xyzi_mod = np.column_stack([in_[:, 0], in_[:, 1],\n",
    "                               in_[:, 2] - z_mod, in_[:, 3]])\n",
    "                               \n",
    "    return xyzi_mod, fitpoints_mod, resid, coeffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# finally, a function to make n-radius smoothed points, using whatever\n",
    "# function works best. Here, mean, median and std are used. Median values\n",
    "# The first two can be used to smooth noisy points. The third is a proxy for\n",
    "# roughness in the neighbourhood chosen.\n",
    "\n",
    "def n_filter(pointcloud, tree, radius):\n",
    "    '''\n",
    "    takes in a point cloud (xyzi), a kDtree (generated in build_tree),\n",
    "    and a neighbourhood.\n",
    "    returns the standard deviation of points within (radius) metres\n",
    "    of each point as a new point list.\n",
    "    '''\n",
    "    nhoods = tree.query_radius(pointcloud[:,0:3], r=radius)\n",
    "    n_stats = []\n",
    "    i = 0\n",
    "    for nhood in nhoods:\n",
    "        #print(nhood)\n",
    "        n_stats.append([np.mean(pointcloud[nhood[:],2]),\\\n",
    "                      np.median(pointcloud[nhood[:],2]),\\\n",
    "                      np.std(pointcloud[nhood[:],2])])\n",
    "        #print(pointcloud[i,:])\n",
    "        #print(new_z[i])\n",
    "        i += 1\n",
    "    return n_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying it all! Setting up reference points and levelling the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### getting some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#setting up a file path for IO\n",
    "infile_dir = '/media/adam/data/is6_f11/lidar.matlab/swaths/'\n",
    "infile = 'is6_f11_pass1_aa_nr2_522816_523019_c.xyz'\n",
    "outfilename = 'is6_f11_pass1_aa_nr2_522816_523019_zi_data.xyz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "input_points = input_points[input_points[:,3] < -5,:]\n",
    "\n",
    "# crop any that are likely to be air returns\n",
    "input_points = input_points[input_points[:,3] < -5,:]\n",
    "\n",
    "input_points = input_points[1:200000,:]\n",
    "\n",
    "#grab an xyzi subset to work with\n",
    "xyzi_ = input_points[:, 1:5]\n",
    "\n",
    "print('x limits: {} - {}'.format(np.min(xyzi_[:,0]), np.max(xyzi_[:,0])))\n",
    "print('y limits: {} - {}'.format(np.min(xyzi_[:,1]), np.max(xyzi_[:,1])))\n",
    "print('total points: {}'.format(len(xyzi_[:,1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "...now we have a set of LiDAR points in local coordinates - oriented such that the X-Y plane is aligned with 'across track', and 'along track' respectively. Essentially pretending the aircraft flew approximately along the Y axis - for interpolation/extrapolation reasons.\n",
    "\n",
    "### next, set up some snow and ice model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#density parameters, should add as a data dictionary\n",
    "d_snow = 326.31 #SIPEX2 snow density (kg/m^3)\n",
    "#d_snow = 305.67 #mean of all EA obs\n",
    "sd_dsnow = 10 #uncertainty in snow density\n",
    "d_ice = 919.6 #empirically derived from matching with AUV draft - kg/m^3, unrealistically high!\n",
    "sd_dice = 10 #uncertainty in ice density\n",
    "d_water = 1028 #Hutchings2015 - seawater density, kg/m^3\n",
    "sd_dwater = 1 #uncertainty in seawater density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#using a loop, recursively iterate fit points toward 'level'\n",
    "\n",
    "i = 0\n",
    "while i <= 2:\n",
    "    #xyzi_, keypoints = fit_points(xyzi_, [10, 2.0], knnf)\n",
    "    xyzi_, keypoints, redi, coeff = f_points(xyzi_, [20, 2], [1 ,5], 10 )\n",
    "    print('keypoints: {}'.format(len(keypoints)))\n",
    "    print('keypoints mean elev: {}'.format(np.mean(keypoints[:,2])))\n",
    "    print('mean elev: {}'.format(np.mean(xyzi_[:,2])))\n",
    "    \n",
    "    print(i)\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##assuming elev > 3m = iceberg!\n",
    "# or 5m here...\n",
    "nobergs = np.where(xyzi_[:,2] < 5)\n",
    "xyzi_ = xyzi_[nobergs[0],:]\n",
    "input_points = input_points[nobergs[0],:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# because water is usually still too high, or too low...\n",
    "# here we collect the lowest values from the fit point set and use them to\n",
    "# adjust point heights. *NB this is *after* the fitting-to-surface process,\n",
    "# and really only required because LiDAR scatters in both directions from\n",
    "# the intended plane of observation.\n",
    "adjust = np.percentile(keypoints[:,2],1)\n",
    "print('adjustment for heights after model fitting: {}'.format(adjust))\n",
    "xyzi_[:,2] = xyzi_[:,2] + adjust"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Build a KDtree. Why? so we can query it for point neighbours! The tree is only built\n",
    "# once, then recycled for as many uses as required. It's a pretty long process\n",
    "print('building KDtree')\n",
    "startTime = datetime.now()\n",
    "points_kdtree = build_tree(xyzi_[:,0:3], 60)\n",
    "print('time to build tree: {}'.format(datetime.now() - startTime))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Query the tree for points within 1 m of each point, and generate mean/median/SD. This\n",
    "# approximates a median smoothing method (or mean smoothing, or...) in a 2m diameter circle\n",
    "# placed around each point\n",
    "startTime = datetime.now()\n",
    "nhood_stats = np.array(n_filter(xyzi_, points_kdtree, 1))\n",
    "print('time to generate n stats: {}'.format(datetime.now() - startTime))\n",
    "\n",
    "#replace xyzi Z with median Z (nhood_stats[:,1])\n",
    "# or mean_z (nhood_stats[:,0])\n",
    "xyzi2 = np.column_stack([input_points[:,1], input_points[:,2], \\\n",
    "                            nhood_stats[:,1], input_points[:,4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Because the tree is already built, collect another dataset smoothed over a bigger radius.\n",
    "# in this case 5.5 m, to make 11 m radius windows at which resolution correlations between \n",
    "# total freeboard and draft begin to make sense (according to Doble et al (2011))\n",
    "startTime = datetime.now()\n",
    "r_proxy = np.array(n_filter(xyzi_, points_kdtree, 5.5 ))\n",
    "print('time to generate n stats: {}'.format(datetime.now() - startTime))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: making sea ice thickness!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#smooth things just a little - assign total freeboard from the 2m diameter nhood filter - in this case median.\n",
    "tf = nhood_stats[:,1]\n",
    "#or keep them raw\n",
    "#tf = xyzi_[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Build traps for points which make no sense.\n",
    "sub_z = np.where(tf<0)\n",
    "print('total points where total freeboard < 0 (reassigned to 0): {}'.format(len(sub_z[0])))\n",
    "tf[sub_z] = 0\n",
    "\n",
    "#assign the lidar point Z uncertainty to total freeboard\n",
    "sd_tf = input_points[:, 8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Estimate some snow\n",
    "zs, zs_uncert = compute_zs(tf, s_i, sd_tf)\n",
    "\n",
    "#...and filter points which make no sense - negative snow values\n",
    "sub_z = np.where(zs<0)\n",
    "print('total points where snow < 0 (reasigned to 0): {}'.format(len(sub_z[0])))\n",
    "zs[sub_z] = 0\n",
    "\n",
    "#...and values where \n",
    "oversnowed = np.where(zs > tf)\n",
    "print('total oversnowed points (zs > tf, zs reassigned to 0): {}'.format(len(oversnowed[0])))\n",
    "zs[oversnowed] = 0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## To do:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- maybe make a class like:\n",
    "    \n",
    "    - pc = pointcloud()\n",
    "    - pc.medianfilter(args)\n",
    "    - pc.tolas(args)\n",
    "    - pc.tonetcdf(args)\n",
    "    - pc.anglefilter(args)\n",
    "    - pc.rangegate(args)\n",
    "    - pc.density(args)\n",
    "    - pc.grid(args)\n",
    "    - pc.normals()\n",
    "    - pc.surf() (maybe)\n",
    "    - pc.outliers()\n",
    "    \n",
    "- and some sea ice specific ones:\n",
    "    - pc.estimatezs(modelparams)\n",
    "    - pc.estimatezi(modelparams)\n",
    "    \n",
    "- and really later:\n",
    "    - pc.make(range, angle, trajectory, transforms)\n",
    "    - pc.computeUncert(args)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}