Merge pull request #2 from jhkennedy/missing-interp

Use a priority interpolation method.

build_greenland_datasets.py

@@ -126,8 +126,9 @@
 #=================================
 speak.notquiet(args,"\nGetting thk, topg, and topgerr from the mass conserving bed data.")
 
-icebridge.get_mcb(args, nc_massCon, nc_base, trans, proj_eigen_gl04c, proj_epsg3413)
+icebridge.get_mcb(args, nc_massCon, nc_bamber, nc_base, base, trans, proj_eigen_gl04c, proj_epsg3413)
 
+nc_bamber.close()
 nc_massCon.close()
 #==== Zurich mask ====
 # apply mask, and get 

data/icebridge.py

@@ -1,67 +1,162 @@
 import pyproj
+import scipy
 import numpy as np
-from scipy import interpolate
 
-from util.ncfunc import copy_atts
 from util import speak
+from util.ncfunc import copy_atts, copy_atts_bad_fill
+from util.projections import DataGrid
+import util.interpolate as interp
 
-def get_mcb(args, nc_massCon, nc_base, trans, proj_eigen_gl04c, proj_epsg3413):
+def get_mcb(args, nc_massCon, nc_bamber, nc_base, base, trans, proj_eigen_gl04c, proj_epsg3413):
     """Get the mass conserving bed data.
     """
-    massCon_y = nc_massCon.variables['y']
-    massCon_ny = massCon_y[:].shape[0]
-
-    massCon_x = nc_massCon.variables['x']
-    massCon_nx = massCon_x[:].shape[0]
-
-    massCon_data = np.ndarray( (massCon_ny,massCon_nx) )
-    base_data    = np.ndarray( (trans.ny,trans.nx) )
-    temp_data    = np.ndarray( (trans.ny,trans.nx) )
-
-    var_list       = [ 'thickness', 'bed',  'errbed' ]
-    rename_massCon = [ 'thk',       'topg', 'topgerr']
-    for vv in range(0, len(var_list) ) :
-        var = var_list[vv]
-        massCon_data[:,:] = 0.
-        base_data[:,:]    = 0.
-        temp_data[:,:]    = 0.
-
-        massCon_var = nc_massCon.variables[var]
-        massCon_data[:,:] = massCon_var[::-1,:] # y fliped when compaired to Bamber
-
-        speak.verbose(args,"   Interpolating "+var_list[vv]+".")
-        massCon_to_base = interpolate.RectBivariateSpline( massCon_y[::-1], massCon_x[:], massCon_data, kx=1, ky=1, s=0) # regular 2d linear interp. but faster
-
-        for ii in range(0, trans.nx):
-            base_data[:,ii] = massCon_to_base.ev(trans.y_grid[:,ii], trans.x_grid[:,ii] )
-
-        # This is only needed for data that is actually referenced from the EIGEN-GL04C Geoid -- all topographical 'z' data. Not needed if not 'z' data.
-        if (var_list[vv] != 'errbed') :
-            temp_x_grid, temp_y_grid, temp_data = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, trans.x_grid.flatten(), trans.y_grid.flatten(), base_data.flatten())
-            temp_data = temp_data.reshape((trans.ny,trans.nx))
-            base_data[:,:] = temp_data[:,:]
-
-        if (var_list[vv] == 'thickness') :
-            bad_val = 0
-        else :
-            bad_val = -9999
-
-        data_min = np.amin( massCon_data[ massCon_data[:,:] != bad_val] )
-        data_max = np.amax( massCon_data[ massCon_data[:,:] != bad_val] )
-
-        data_mask = np.ma.masked_outside(base_data, data_min, data_max)
-        base_data[ data_mask.mask ] = bad_val
-
-        speak.verbose(args,"   Writing "+rename_massCon[vv]+" to base.")
-        base_var = nc_base.createVariable( rename_massCon[vv], 'f4', ('y','x',) )
-        base_var[:,:] = base_data[:,:]
-        #copy_atts(massCon_var, base_var)
-        #NOTE: this data is formated short, while the rest of the data is a float. Does not like _FillValue from short data. 
-        atts = massCon_var.ncattrs()
-        for ii in range(len(atts)):
-            if (atts[ii] != '_FillValue'):
-                base_var.setncattr(atts[ii], massCon_var.getncattr(atts[ii]))
+    massCon = DataGrid()
+
+    massCon.y = nc_massCon.variables['y']
+    massCon.x = nc_massCon.variables['x']
+    massCon.ny = massCon.y[:].shape[0]
+    massCon.nx = massCon.x[:].shape[0]
+    massCon.make_grid_flip_y()
+
+    massCon.yx = np.ndarray( (len(massCon.y_grid.ravel()),2) )
+    massCon.yx[:,0] = massCon.y_grid.ravel()
+    massCon.yx[:,1] = massCon.x_grid.ravel()
+
+    massCon.tree = scipy.spatial.cKDTree(massCon.yx)
+
+    trans.yx = np.ndarray( (len(trans.y_grid.ravel()),2) )
+    trans.yx[:,0] = trans.y_grid.ravel()
+    trans.yx[:,1] = trans.x_grid.ravel()
+
+    trans.qd, trans.qi = massCon.tree.query(trans.yx, k=1)   # nearest neighbor in massCon for transformed base grid
+
+    massCon.thickness = nc_massCon.variables['thickness']
+    massCon.thk = np.ndarray(massCon.dims)
+    massCon.thk[:,:] = massCon.thickness[::-1,:] # y fliped when compaired to Bamber
+
+    speak.verbose(args,"   Interpolating thickness.")
+    massCon_to_base = scipy.interpolate.RectBivariateSpline( massCon.y[::-1], massCon.x[:], massCon.thk, kx=1, ky=1, s=0) # regular 2d linear interp. but faster
+    trans.thk = np.zeros( trans.dims )
+    for ii in range(0, base.nx):
+        trans.thk[:,ii] = massCon_to_base.ev(trans.y_grid[:,ii], trans.x_grid[:,ii] )
+
+
+    speak.verbose(args,"   Writing thk to base.")
+    base.thk = nc_base.createVariable('thk', 'f4', ('y','x',) )
+    base.thk[:,:] = trans.thk[:,:]
+    copy_atts_bad_fill(massCon.thickness, base.thk, -9999.)
+
+    speak.verbose(args,"   Interpolating, with priority, topg and topgerr.")
+    speak.verbose(args,"      Primary Data [IceBridge]:  bed and errbed.")
+    speak.verbose(args,"      Secondary Data [BamberDEM]:  BedrockElevation and BedrockError.")
+
+    pri_data = np.ma.masked_equal( nc_massCon.variables['bed'][::-1,:] , -9999)
+    sec_data = np.ma.masked_values( nc_bamber.variables['BedrockElevation'][:,:], -9999.)
+    new_data = np.ma.array(np.zeros(base.dims), mask=np.zeros(base.dims))
+
+    pri_err = np.ma.masked_equal( nc_massCon.variables['errbed'][::-1,:] , -9999)
+    sec_err = np.ma.masked_values( nc_bamber.variables['BedrockError'][:,:], -9999.)
+    new_err = np.ma.array(np.zeros(base.dims), mask=np.zeros(base.dims))
+
+    for ii in range(0, base.ny):
+        for jj in range(0, base.nx):
+            # make sure inside priority grid (massCon)
+            if (trans.y_grid[ii,jj] < massCon.y_grid[0,0] or trans.y_grid[ii,jj] > massCon.y_grid[-1,0]):
+                # outside y range
+                if sec_data.mask[ii,jj]:
+                    new_data.mask[ii,jj] = True
+                    new_err.mask[ii,jj] = True
+                    continue
+                else:
+                    tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, base.x_grid[ii,jj], base.y_grid[ii,jj], sec_data[ii,jj])
+                    new_data[ii,jj] = td
+                    new_err[ii,jj] = sec_err[ii,jj]
+                    continue
+
+            if (trans.x_grid[ii,jj] < massCon.x_grid[0,0] or trans.x_grid[ii,jj] > massCon.x_grid[0,-1]):
+                # outside x range 
+                if sec_data.mask[ii,jj]:
+                    new_data.mask[ii,jj] = True
+                    new_err.mask[ii,jj] = True
+                    continue
+                else:
+                    tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, base.x_grid[ii,jj], base.y_grid[ii,jj], sec_data[ii,jj])
+                    new_data[ii,jj] = td
+                    new_err[ii,jj] = sec_err[ii,jj]
+                    continue
+
+            indx = np.ravel_multi_index( (ii,jj), base.dims)
+
+            # nearest neighbor indices
+            nn_ii, nn_jj = np.unravel_index( trans.qi[indx], massCon.dims)
+
+            # to find nearest neighbors quadrent
+            i_s = -1
+            j_s = -1
+
+            # find quadrent
+            if trans.y_grid[ii,jj] >= massCon.y_grid[nn_ii,nn_jj]:
+                i_s = +1
+            if trans.x_grid[ii,jj] >= massCon.x_grid[nn_ii,nn_jj]:
+                j_s = +1
+
+            # check for missing priority data!
+            missing_points, interp_dict = interp.check_missing(pri_data, (nn_ii, nn_jj), i_s, j_s)
+            missing_err_pts, err_dict = interp.check_missing(pri_err, (nn_ii, nn_jj), i_s, j_s)
+
+            # get secondary data!
+            if not missing_points :
+                pass
+
+            elif not sec_data.mask[ii,jj] :
+                tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, base.x_grid[ii,jj], base.y_grid[ii,jj], sec_data[ii,jj])
+                if len(missing_points) <= 3 :
+                    for point in missing_points:
+                        # use secondary data at (ii,jj) for missing points, but keep same interp weight!
+                        interp_dict[point] = td
+                        err_dict[point] = sec_err[ii,jj]
+
+                else:
+                    new_data[ii,jj] = td
+                    new_err[ii,jj] = sec_err[ii,jj]
+                    continue
+
             else:
-                base_var.setncattr('missing_value', base_data[-1,-1]) # from a known bad point
+                new_data.mask[ii,jj] = True
+                new_err.mask[ii,jj] = True
+                continue
+
+            # interpolate!
+            alpha = ( trans.y_grid[ii,jj] - massCon.y_grid[nn_ii,nn_jj] )/(massCon.dy*i_s)
+            beta  = ( trans.x_grid[ii,jj] - massCon.x_grid[nn_ii,nn_jj] )/(massCon.dx*j_s)
+
+            w = interp.linear_weights(alpha, beta)
+
+            new_data[ii,jj] = interp_dict[ (nn_ii,    nn_jj    ) ]*w[0] \
+                             +interp_dict[ (nn_ii,    nn_jj+j_s) ]*w[1] \
+                             +interp_dict[ (nn_ii+i_s,nn_jj+j_s) ]*w[2] \
+                             +interp_dict[ (nn_ii+i_s,nn_jj    ) ]*w[3]
+
+
+            new_err[ii,jj] = err_dict[ (nn_ii,    nn_jj    ) ]*w[0] \
+                            +err_dict[ (nn_ii,    nn_jj+j_s) ]*w[1] \
+                            +err_dict[ (nn_ii+i_s,nn_jj+j_s) ]*w[2] \
+                            +err_dict[ (nn_ii+i_s,nn_jj    ) ]*w[3]
+
+            missing_points = None
+            interp_dict = None
+            err_dict = None
+    # now transform new data back to bamber grid. 
+    temp_x_grid, temp_y_grid, temp_data = pyproj.transform(proj_epsg3413, proj_eigen_gl04c, trans.x_grid.flatten(), trans.y_grid.flatten(), new_data.flatten())
+    new_data[:,:] = temp_data.reshape( base.dims )[:,:]
+
+    speak.verbose(args,"   Writing topg topgerr to base.")
+    base.topg = nc_base.createVariable('topg', 'f4', ('y','x',) )
+    base.topg[:,:] = new_data.filled(-9999.)[:,:]
+    copy_atts_bad_fill(nc_massCon.variables['bed'], base.topg, -9999.)
 
+    base.topgerr = nc_base.createVariable('topgerr', 'f4', ('y','x',) )
+    base.topgerr[:,:] = new_err.filled(-9999.)[:,:]
+    copy_atts_bad_fill(nc_massCon.variables['errbed'], base.topgerr, -9999.)
+
 

play_tree.py

@@ -0,0 +1,29 @@
+import numpy as np
+import scipy
+from scipy.spatial import cKDTree
+
+x = [1,2,3]
+y = [1,2,3]
+
+xx, yy =  scipy.meshgrid(x, y, indexing='ij')
+
+xy = np.ndarray( (len(xx.ravel() ), 2) )
+
+xy[:,0] = xx.ravel()
+xy[:,1] = yy.ravel()
+
+print(xy)
+
+tree = cKDTree( xy )
+
+qp = (4., 4.)
+d, indx = tree.query( qp, k=4) 
+
+print('x neighbors: ')
+print(xy[indx,0])
+
+print('y neighbors: ')
+print(xy[indx,1])
+
+print('Distances: ')
+print(d)

util/interpolate.py

@@ -0,0 +1,32 @@
+import numpy as np
+import scipy
+
+def check_missing(data, anchor, i_s, j_s):
+    missing_points = []
+    interp_dict = {}
+
+    dims_list = [(anchor[0],    anchor[1]    ),
+                 (anchor[0],    anchor[1]+j_s),
+                 (anchor[0]+i_s,anchor[1]+j_s),
+                 (anchor[0]+i_s,anchor[1]    )]
+
+    for dims in dims_list :
+        if not any( d < 0 for d in dims) and data.mask[dims]:
+            missing_points.append(dims)
+
+        interp_dict.update({dims: data[dims]} )
+
+
+    return (missing_points, interp_dict)
+
+
+def linear_weights(alpha, beta):
+
+    w = np.zeros(4)
+
+    w[0] = (1.-alpha)*(1.-beta)
+    w[1] = (1.-alpha)*(   beta)
+    w[2] = (   alpha)*(   beta)
+    w[3] = (   alpha)*(1.-beta)
+
+    return w

util/ncfunc.py

@@ -46,3 +46,37 @@ def copy_atts(fin, fout) :
         fout.setncattr(atts[ii], fin.getncattr(atts[ii]))
 
 
+def copy_atts_bad_fill(fin, fout, missing_value) :
+    """Copy all netCDF attributes except _FillValue.  
+
+    This function copies all the attributes from one netCDF element to another,
+    but ignores the _FillValue attribute and sets MissingValue. 
+    
+    Parameters
+    ----------
+    fin : 
+        Source netCDF element
+    fout :
+        Target netCDF element
+    missing_value :
+        Value to set as indicator of missing values.
+
+    Examples
+    --------
+    Copy the attributes from one variable to another.
+
+    >>> old_var = nc_old.variables['old']
+    >>> new_var = nc_new.createVariable('new', 'f4', ('y','x',) )
+    >>> new_var[:,:] = old_var[:,:]
+    >>> copy_atts_bad_fill( old_var,new_var, -9999. )
+    """
+
+    # get a list of global attribute names from the incoming file
+    atts = fin.ncattrs()
+
+    # place those attributes in the outgoing file
+    for ii in range(len(atts)) :
+        if (atts[ii] != '_FillValue'):
+            fout.setncattr(atts[ii], fin.getncattr(atts[ii]))
+        else:
+            fout.setncattr('missing_value', missing_value)

util/projections.py

@@ -12,6 +12,17 @@ def make_grid(self):
         """A way to make a basic grid from x,y data.
         """
         self.y_grid, self.x_grid = scipy.meshgrid(self.y[:], self.x[:], indexing='ij')
+        self.dims = (self.ny, self.nx)
+        self.dy = self.y[1]-self.y[0]
+        self.dx = self.x[1]-self.x[0]
+
+    def make_grid_flip_y(self):
+        """A way to make a basic grid from x,y data, inverting y.
+        """
+        self.y_grid, self.x_grid = scipy.meshgrid(self.y[::-1], self.x[:], indexing='ij')
+        self.dims = (self.ny, self.nx)
+        self.dy = self.y[0]-self.y[1]
+        self.dx = self.x[1]-self.x[0]
 
 
 def greenland(args, lc_bamber, base):
@@ -40,6 +51,7 @@ def greenland(args, lc_bamber, base):
     trans = DataGrid()
     trans.ny = base.ny
     trans.nx = base.nx
+    trans.dims = base.dims
 
     trans.x_grid, trans.y_grid = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, base.x_grid.flatten(), base.y_grid.flatten())
     trans.y_grid = trans.y_grid.reshape((base.ny,base.nx))