From 0c9cb750c007f2cd04997df57cdae2e1f519be6a Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Thu, 5 Mar 2015 18:21:50 -0500
Subject: [PATCH 1/7] Start playing with kdtree

---
 test_tree.py        | 30 ++++++++++++++++++++++++++++++
 util/interpolate.py |  4 ++++
 2 files changed, 34 insertions(+)
 create mode 100644 test_tree.py
 create mode 100644 util/interpolate.py

diff --git a/test_tree.py b/test_tree.py
new file mode 100644
index 0000000..331d686
--- /dev/null
+++ b/test_tree.py
@@ -0,0 +1,30 @@
+
+import numpy as np
+import scipy
+
+from util.ncfunc import get_nc_file
+from util.projections import DataGrid
+
+lc_bamber   = 'data/BamberDEM/Greenland_bedrock_topography_V3.nc'
+lc_massCon  = 'data/IceBridge/Greenland/MCdataset-2014-11-19.nc'
+
+
+nc_bamber = get_nc_file(lc_bamber,'r')
+nc_massCon = get_nc_file(lc_massCon,'r')
+
+bamber  = DataGrid()
+massCon = DataGrid()
+
+bamber.y = nc_bamber.variables['projection_y_coordinate']
+bamber.x = nc_bamber.variables['projection_x_coordinate']
+bamber.ny = bamber.y[:].shape[0]
+bamber.nx = bamber.x[:].shape[0]
+bamber.make_grid()
+
+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()
+
+bamber.r_grid = ( bamber.y_grid.ravel(), bamber.x_grid.ravel() )
diff --git a/util/interpolate.py b/util/interpolate.py
new file mode 100644
index 0000000..4bd5dbf
--- /dev/null
+++ b/util/interpolate.py
@@ -0,0 +1,4 @@
+import numpy as np
+import scipy
+
+

From 436c92a5902c5f41d018ea7f0722f622e2993791 Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Mon, 9 Mar 2015 18:30:55 -0400
Subject: [PATCH 2/7] Began building interpolation method

---
 test_tree.py | 56 +++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 53 insertions(+), 3 deletions(-)

diff --git a/test_tree.py b/test_tree.py
index 331d686..2e2ef8c 100644
--- a/test_tree.py
+++ b/test_tree.py
@@ -1,6 +1,8 @@
-
-import numpy as np
+import os
 import scipy
+import pyproj
+import numpy as np
+from scipy.spatial import cKDTree
 
 from util.ncfunc import get_nc_file
 from util.projections import DataGrid
@@ -8,6 +10,7 @@
 lc_bamber   = 'data/BamberDEM/Greenland_bedrock_topography_V3.nc'
 lc_massCon  = 'data/IceBridge/Greenland/MCdataset-2014-11-19.nc'
 
+path_bamber = os.path.dirname(lc_bamber)
 
 nc_bamber = get_nc_file(lc_bamber,'r')
 nc_massCon = get_nc_file(lc_massCon,'r')
@@ -21,10 +24,57 @@
 bamber.nx = bamber.x[:].shape[0]
 bamber.make_grid()
 
+bamber.yx = np.ndarray( (len( bamber.y_grid.ravel() ),2) )
+bamber.yx[:,0] = bamber.y_grid.ravel()
+bamber.yx[:,1] = bamber.x_grid.ravel()
+
 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()
 
-bamber.r_grid = ( bamber.y_grid.ravel(), bamber.x_grid.ravel() )
+massCon.yx = np.ndarray( (len( massCon.y_grid.ravel() ),2) )
+massCon.yx[:,0] = massCon.y_grid.ravel()
+massCon.yx[:,1] = massCon.x_grid.ravel()
+
+# ckd trees of data grids
+bamber.tree = cKDTree(bamber.yx)
+massCon.tree = cKDTree(massCon.yx)
+
+# get transformation grids
+proj_eigen_gl04c = pyproj.Proj('+proj=stere +lat_ts=71.0 +lat_0=90 +lon_0=321.0 +k_0=1.0 +geoidgrids='+path_bamber+'/egm08_25.gtx')
+proj_epsg3413 = pyproj.Proj('+proj=stere +lat_ts=70.0 +lat_0=90 +lon_0=-45.0 +k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84 +units=m')
+
+b_trans = DataGrid()
+b_trans.ny = bamber.ny
+b_trans.nx = bamber.nx
+
+b_trans.x_grid, b_trans.y_grid = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid.flatten(), bamber.y_grid.flatten())
+b_trans.y_grid = b_trans.y_grid.reshape((b_trans.ny,b_trans.nx))
+b_trans.x_grid = b_trans.x_grid.reshape((b_trans.ny,b_trans.nx))
+
+b_trans.yx = np.ndarray( (len( b_trans.y_grid.ravel() ),2) )
+b_trans.yx[:,0] = b_trans.y_grid.ravel()
+b_trans.yx[:,1] = b_trans.x_grid.ravel()
+
+m_trans = DataGrid()
+m_trans.ny = massCon.ny
+m_trans.nx = massCon.nx
+
+m_trans.x_grid, m_trans.y_grid = pyproj.transform(proj_epsg3413, proj_eigen_gl04c, massCon.x_grid.flatten(), massCon.y_grid.flatten())
+m_trans.y_grid = m_trans.y_grid.reshape((m_trans.ny,m_trans.nx))
+m_trans.x_grid = m_trans.x_grid.reshape((m_trans.ny,m_trans.nx))
+
+m_trans.yx = np.ndarray( (len( m_trans.y_grid.ravel() ),2) )
+m_trans.yx[:,0] = m_trans.y_grid.ravel()
+m_trans.yx[:,1] = m_trans.x_grid.ravel()
+
+# find transform grids nearest neighbor
+b_trans.qd, b_trans.qi = massCon.tree.query(b_trans.yx, k=1)
+m_trans.qd, m_trans.qi = bamber.tree.query(m_trans.yx, k=1)
+
+# build data array.
+new_data = np.ma.array( np.zeros( (bamber.ny,bamber.nx) ), mask=np.zeros( (bamber.ny,bamber.nx) ) )
+
+

From 8ad7900b16b3cfcae02eb03a774702f4202ccf46 Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Mon, 9 Mar 2015 18:31:25 -0400
Subject: [PATCH 3/7] A script for playing with cKDTrees.

---
 play_tree.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 play_tree.py

diff --git a/play_tree.py b/play_tree.py
new file mode 100644
index 0000000..0b0ffb7
--- /dev/null
+++ b/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)

From f25f364f109c497b50778e74242b9847f37cf4f2 Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Wed, 11 Mar 2015 15:17:57 -0400
Subject: [PATCH 4/7] Finished interpolation method that uses primary and
 secondary data

Will first try to use the primary data for linear interpolation.

If there are missing values in the primary set, it will use the
secondary data.
---
 test_tree.py        | 148 ++++++++++++++++++++++++++++++++++++++------
 util/interpolate.py |  28 +++++++++
 util/projections.py |  11 ++++
 3 files changed, 167 insertions(+), 20 deletions(-)

diff --git a/test_tree.py b/test_tree.py
index 2e2ef8c..a741b69 100644
--- a/test_tree.py
+++ b/test_tree.py
@@ -2,11 +2,16 @@
 import scipy
 import pyproj
 import numpy as np
+from netCDF4 import Dataset
 from scipy.spatial import cKDTree
 
-from util.ncfunc import get_nc_file
+import util.interpolate as interp
+from util.ncfunc import copy_atts, get_nc_file
 from util.projections import DataGrid
 
+
+print('Loading Data')
+
 lc_bamber   = 'data/BamberDEM/Greenland_bedrock_topography_V3.nc'
 lc_massCon  = 'data/IceBridge/Greenland/MCdataset-2014-11-19.nc'
 
@@ -32,16 +37,34 @@
 massCon.x = nc_massCon.variables['x']
 massCon.ny = massCon.y[:].shape[0]
 massCon.nx = massCon.x[:].shape[0]
-massCon.make_grid()
+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()
 
+print('Creating test data file')
+# create some base variables
+class test():
+    pass
+
+nc_test = Dataset('test.nc','w')
+nc_test.createDimension('y',bamber.ny)
+nc_test.createDimension('x',bamber.nx)
+test.y = nc_test.createVariable('y', 'f4', 'y')
+test.y[:] = bamber.y[:]
+copy_atts(bamber.y, test.y) #FIXME: units say km, but it's actuall in m
+
+test.x = nc_test.createVariable('x', 'f4', 'x')
+test.x[:] = bamber.x[:]
+copy_atts(bamber.x, test.x) #FIXME: units say km, but it's actuall in m
+
+print('Building trees')
 # ckd trees of data grids
 bamber.tree = cKDTree(bamber.yx)
 massCon.tree = cKDTree(massCon.yx)
 
+print('Transformations')
 # get transformation grids
 proj_eigen_gl04c = pyproj.Proj('+proj=stere +lat_ts=71.0 +lat_0=90 +lon_0=321.0 +k_0=1.0 +geoidgrids='+path_bamber+'/egm08_25.gtx')
 proj_epsg3413 = pyproj.Proj('+proj=stere +lat_ts=70.0 +lat_0=90 +lon_0=-45.0 +k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84 +units=m')
@@ -49,6 +72,7 @@
 b_trans = DataGrid()
 b_trans.ny = bamber.ny
 b_trans.nx = bamber.nx
+b_trans.dims = (bamber.ny, bamber.nx)
 
 b_trans.x_grid, b_trans.y_grid = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid.flatten(), bamber.y_grid.flatten())
 b_trans.y_grid = b_trans.y_grid.reshape((b_trans.ny,b_trans.nx))
@@ -58,23 +82,107 @@
 b_trans.yx[:,0] = b_trans.y_grid.ravel()
 b_trans.yx[:,1] = b_trans.x_grid.ravel()
 
-m_trans = DataGrid()
-m_trans.ny = massCon.ny
-m_trans.nx = massCon.nx
+b_trans.qd, b_trans.qi = massCon.tree.query(b_trans.yx, k=1)   # nearest neighbor in massCon for transformed bamber grid
 
-m_trans.x_grid, m_trans.y_grid = pyproj.transform(proj_epsg3413, proj_eigen_gl04c, massCon.x_grid.flatten(), massCon.y_grid.flatten())
-m_trans.y_grid = m_trans.y_grid.reshape((m_trans.ny,m_trans.nx))
-m_trans.x_grid = m_trans.x_grid.reshape((m_trans.ny,m_trans.nx))
-
-m_trans.yx = np.ndarray( (len( m_trans.y_grid.ravel() ),2) )
-m_trans.yx[:,0] = m_trans.y_grid.ravel()
-m_trans.yx[:,1] = m_trans.x_grid.ravel()
-
-# find transform grids nearest neighbor
-b_trans.qd, b_trans.qi = massCon.tree.query(b_trans.yx, k=1)
-m_trans.qd, m_trans.qi = bamber.tree.query(m_trans.yx, k=1)
-
-# build data array.
+print('Bulding data arrays')
+# build data arrays
+pri_data = np.ma.masked_equal( nc_massCon.variables['bed'][::-1,:], -9999)     # MASSCON GRID
+sec_data = np.ma.masked_values( nc_bamber.variables['BedrockElevation'][:,:], -9999.) # BAMBER GRID
 new_data = np.ma.array( np.zeros( (bamber.ny,bamber.nx) ), mask=np.zeros( (bamber.ny,bamber.nx) ) )
-
-
+reason = np.zeros( bamber.dims )
+
+#import matplotlib.pyplot as plt
+#plt.imshow(pri_data.filled(-1))
+#plt.show()
+
+print('Interpolating...')
+# loop through all data points
+for ii in range(0, bamber.ny):
+    for jj in range(0, bamber.nx):
+        # interpolate bamber point on massCon grid
+        indx = np.ravel_multi_index( (ii,jj),b_trans.dims )
+
+        nn_mc = b_trans.qi[indx]
+        nn_ii, nn_jj = np.unravel_index(nn_mc, massCon.dims)
+        
+        y_s = -1
+        x_s = -1
+        
+        # make sure inside priority grid (massCon)
+        if (b_trans.y_grid[ii,jj] < massCon.y_grid[0,0] or b_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
+                reason[ii,jj] = 1
+                continue
+            else:
+                tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid[ii,jj], bamber.y_grid[ii,jj], sec_data[ii,jj])
+                new_data[ii,jj] = td
+                reason[ii,jj] = -1
+                continue
+        
+        if (b_trans.x_grid[ii,jj] < massCon.x_grid[0,0] or b_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
+                reason[ii,jj] = 2
+                continue
+            else:
+                tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid[ii,jj], bamber.y_grid[ii,jj], sec_data[ii,jj])
+                new_data[ii,jj] = td
+                reason[ii,jj] = -1
+                continue
+
+        # find quadrent
+        if b_trans.y_grid[ii,jj] >= massCon.y_grid[nn_ii,nn_jj]:
+            y_s = +1
+        if b_trans.x_grid[ii,jj] >= massCon.x_grid[nn_ii,nn_jj]:
+            x_s = +1
+
+        # check for missing priority data!
+        missing_points, interp_dict = interp.check_missing(pri_data, (nn_ii, nn_jj), y_s, x_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, bamber.x_grid[ii,jj], bamber.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
+                reason[ii,jj] = -2
+                continue
+            else:
+                new_data[ii,jj] = td
+                reason[ii,jj] = -1
+                continue
+        else:
+            new_data.mask[ii,jj] = True
+            reason[ii,jj] = 3
+            continue
+
+        # interpolate!
+        alpha = ( b_trans.y_grid[ii,jj] - massCon.y_grid[nn_ii,nn_jj] )/(massCon.dy*y_s)
+        beta  = ( b_trans.x_grid[ii,jj] - massCon.x_grid[nn_ii,nn_jj] )/(massCon.dx*x_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+x_s) ]*w[1] \
+                         +interp_dict[ (nn_ii+y_s,nn_jj+x_s) ]*w[2] \
+                         +interp_dict[ (nn_ii+y_s,nn_jj    ) ]*w[3]
+
+import matplotlib.pyplot as plt
+plt.imshow(new_data.mask)
+plt.show()
+
+print('Writing data')
+test_var = nc_test.createVariable( 'interpData', 'f4', ('y','x',) )
+test_var[:,:] = new_data[:,:]
+
+test_reason = nc_test.createVariable('reason','i',('y','x',))
+test_reason[:,:] = reason[:,:]
+nc_test.close()
+print('And Done!')
diff --git a/util/interpolate.py b/util/interpolate.py
index 4bd5dbf..c26ff5a 100644
--- a/util/interpolate.py
+++ b/util/interpolate.py
@@ -1,4 +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
diff --git a/util/projections.py b/util/projections.py
index af945d8..71d4fb2 100644
--- a/util/projections.py
+++ b/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):

From 38f3f2651e59eb423277c2a6c863363ad56eaf78 Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Wed, 11 Mar 2015 17:15:11 -0400
Subject: [PATCH 5/7] Added a tree plot and fixed errant continue statements

---
 test_tree.py |  8 ++++----
 tree_plot.py | 37 +++++++++++++++++++++++++++++++++++++
 2 files changed, 41 insertions(+), 4 deletions(-)
 create mode 100644 tree_plot.py

diff --git a/test_tree.py b/test_tree.py
index a741b69..89b6b34 100644
--- a/test_tree.py
+++ b/test_tree.py
@@ -145,7 +145,7 @@ class test():
         # 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, bamber.x_grid[ii,jj], bamber.y_grid[ii,jj], sec_data[ii,jj])
             if len(missing_points) <= 3 :
@@ -153,11 +153,11 @@ class test():
                     # use secondary data at (ii,jj) for missing points, but keep same interp weight!
                     interp_dict[point] = td
                 reason[ii,jj] = -2
-                continue
+            
             else:
                 new_data[ii,jj] = td
                 reason[ii,jj] = -1
-                continue
+
         else:
             new_data.mask[ii,jj] = True
             reason[ii,jj] = 3
@@ -175,7 +175,7 @@ class test():
                          +interp_dict[ (nn_ii+y_s,nn_jj    ) ]*w[3]
 
 import matplotlib.pyplot as plt
-plt.imshow(new_data.mask)
+plt.imshow(reason)
 plt.show()
 
 print('Writing data')
diff --git a/tree_plot.py b/tree_plot.py
new file mode 100644
index 0000000..0f636e5
--- /dev/null
+++ b/tree_plot.py
@@ -0,0 +1,37 @@
+import numpy as np
+from util.ncfunc import copy_atts, get_nc_file
+from util.projections import DataGrid
+
+
+nc_test = get_nc_file('test.nc','r')
+test = DataGrid()
+
+test.reason = nc_test.variables['reason']
+test.dims = test.reason[:,:].shape
+
+reason = np.ndarray(test.dims)
+reason[:,:] = test.reason[:,:]
+
+test.img = np.ones( (test.dims[0], test.dims[1], 3) )
+test.img *= 255
+
+# use secondary as interp data! 
+test.img[reason == -2, : ] = [255,0,0]
+# use secondary as new data! 
+test.img[reason == -1, : ] = [102,255,102]
+# outside y range and no secondary data
+test.img[reason == 1, : ] = [255,102,255]
+# outside x range and no secondary data 
+test.img[reason == 2, : ] = [178,102,255]
+# inside but no secondary data
+test.img[reason == 3, : ] = [255,102,178]
+
+
+
+test.img /= 255
+
+import matplotlib.pyplot as plt
+plt.imshow(test.img[::-1,:,:])
+plt.show()
+
+nc_test.close

From caa1dd05520e2ceadb3ca6b82d35aa91f5ad97b4 Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Thu, 12 Mar 2015 11:58:18 -0400
Subject: [PATCH 6/7] Build now uses priority interpolation method!

---
 build_greenland_datasets.py |   3 +-
 data/icebridge.py           | 209 ++++++++++++++++++++++++++----------
 test_tree.py                |   8 +-
 tree_plot.py                |   3 +-
 util/ncfunc.py              |  34 ++++++
 util/projections.py         |   1 +
 6 files changed, 195 insertions(+), 63 deletions(-)

diff --git a/build_greenland_datasets.py b/build_greenland_datasets.py
index 2ac9a9c..671741d 100644
--- a/build_greenland_datasets.py
+++ b/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 
diff --git a/data/icebridge.py b/data/icebridge.py
index f78ff68..f220fb4 100644
--- a/data/icebridge.py
+++ b/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.)
+    
 
diff --git a/test_tree.py b/test_tree.py
index 89b6b34..564d5a2 100644
--- a/test_tree.py
+++ b/test_tree.py
@@ -61,7 +61,6 @@ class test():
 
 print('Building trees')
 # ckd trees of data grids
-bamber.tree = cKDTree(bamber.yx)
 massCon.tree = cKDTree(massCon.yx)
 
 print('Transformations')
@@ -152,11 +151,12 @@ class test():
                 for point in missing_points:
                     # use secondary data at (ii,jj) for missing points, but keep same interp weight!
                     interp_dict[point] = td
-                reason[ii,jj] = -2
+                reason[ii,jj] = -3
             
             else:
                 new_data[ii,jj] = td
-                reason[ii,jj] = -1
+                reason[ii,jj] = -2
+                continue
 
         else:
             new_data.mask[ii,jj] = True
@@ -175,7 +175,7 @@ class test():
                          +interp_dict[ (nn_ii+y_s,nn_jj    ) ]*w[3]
 
 import matplotlib.pyplot as plt
-plt.imshow(reason)
+plt.imshow(reason[::-1,:])
 plt.show()
 
 print('Writing data')
diff --git a/tree_plot.py b/tree_plot.py
index 0f636e5..0955a18 100644
--- a/tree_plot.py
+++ b/tree_plot.py
@@ -16,9 +16,10 @@
 test.img *= 255
 
 # use secondary as interp data! 
-test.img[reason == -2, : ] = [255,0,0]
+test.img[reason == -3, : ] = [255,0,0]
 # use secondary as new data! 
 test.img[reason == -1, : ] = [102,255,102]
+test.img[reason == -2, : ] = [102,255,102]
 # outside y range and no secondary data
 test.img[reason == 1, : ] = [255,102,255]
 # outside x range and no secondary data 
diff --git a/util/ncfunc.py b/util/ncfunc.py
index 980dc96..1a6d755 100644
--- a/util/ncfunc.py
+++ b/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)
diff --git a/util/projections.py b/util/projections.py
index 71d4fb2..7248575 100644
--- a/util/projections.py
+++ b/util/projections.py
@@ -51,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))

From 1b1305e2c510ddcae23674f6c69c95b28365da52 Mon Sep 17 00:00:00 2001
From: Joseph H Kennedy <kennedyjh@ornl.gov>
Date: Thu, 12 Mar 2015 12:01:40 -0400
Subject: [PATCH 7/7] Stop tracking unessesary test files.

---
 test_tree.py | 188 ---------------------------------------------------
 tree_plot.py |  38 -----------
 2 files changed, 226 deletions(-)
 delete mode 100644 test_tree.py
 delete mode 100644 tree_plot.py

diff --git a/test_tree.py b/test_tree.py
deleted file mode 100644
index 564d5a2..0000000
--- a/test_tree.py
+++ /dev/null
@@ -1,188 +0,0 @@
-import os
-import scipy
-import pyproj
-import numpy as np
-from netCDF4 import Dataset
-from scipy.spatial import cKDTree
-
-import util.interpolate as interp
-from util.ncfunc import copy_atts, get_nc_file
-from util.projections import DataGrid
-
-
-print('Loading Data')
-
-lc_bamber   = 'data/BamberDEM/Greenland_bedrock_topography_V3.nc'
-lc_massCon  = 'data/IceBridge/Greenland/MCdataset-2014-11-19.nc'
-
-path_bamber = os.path.dirname(lc_bamber)
-
-nc_bamber = get_nc_file(lc_bamber,'r')
-nc_massCon = get_nc_file(lc_massCon,'r')
-
-bamber  = DataGrid()
-massCon = DataGrid()
-
-bamber.y = nc_bamber.variables['projection_y_coordinate']
-bamber.x = nc_bamber.variables['projection_x_coordinate']
-bamber.ny = bamber.y[:].shape[0]
-bamber.nx = bamber.x[:].shape[0]
-bamber.make_grid()
-
-bamber.yx = np.ndarray( (len( bamber.y_grid.ravel() ),2) )
-bamber.yx[:,0] = bamber.y_grid.ravel()
-bamber.yx[:,1] = bamber.x_grid.ravel()
-
-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()
-
-print('Creating test data file')
-# create some base variables
-class test():
-    pass
-
-nc_test = Dataset('test.nc','w')
-nc_test.createDimension('y',bamber.ny)
-nc_test.createDimension('x',bamber.nx)
-test.y = nc_test.createVariable('y', 'f4', 'y')
-test.y[:] = bamber.y[:]
-copy_atts(bamber.y, test.y) #FIXME: units say km, but it's actuall in m
-
-test.x = nc_test.createVariable('x', 'f4', 'x')
-test.x[:] = bamber.x[:]
-copy_atts(bamber.x, test.x) #FIXME: units say km, but it's actuall in m
-
-print('Building trees')
-# ckd trees of data grids
-massCon.tree = cKDTree(massCon.yx)
-
-print('Transformations')
-# get transformation grids
-proj_eigen_gl04c = pyproj.Proj('+proj=stere +lat_ts=71.0 +lat_0=90 +lon_0=321.0 +k_0=1.0 +geoidgrids='+path_bamber+'/egm08_25.gtx')
-proj_epsg3413 = pyproj.Proj('+proj=stere +lat_ts=70.0 +lat_0=90 +lon_0=-45.0 +k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84 +units=m')
-
-b_trans = DataGrid()
-b_trans.ny = bamber.ny
-b_trans.nx = bamber.nx
-b_trans.dims = (bamber.ny, bamber.nx)
-
-b_trans.x_grid, b_trans.y_grid = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid.flatten(), bamber.y_grid.flatten())
-b_trans.y_grid = b_trans.y_grid.reshape((b_trans.ny,b_trans.nx))
-b_trans.x_grid = b_trans.x_grid.reshape((b_trans.ny,b_trans.nx))
-
-b_trans.yx = np.ndarray( (len( b_trans.y_grid.ravel() ),2) )
-b_trans.yx[:,0] = b_trans.y_grid.ravel()
-b_trans.yx[:,1] = b_trans.x_grid.ravel()
-
-b_trans.qd, b_trans.qi = massCon.tree.query(b_trans.yx, k=1)   # nearest neighbor in massCon for transformed bamber grid
-
-print('Bulding data arrays')
-# build data arrays
-pri_data = np.ma.masked_equal( nc_massCon.variables['bed'][::-1,:], -9999)     # MASSCON GRID
-sec_data = np.ma.masked_values( nc_bamber.variables['BedrockElevation'][:,:], -9999.) # BAMBER GRID
-new_data = np.ma.array( np.zeros( (bamber.ny,bamber.nx) ), mask=np.zeros( (bamber.ny,bamber.nx) ) )
-reason = np.zeros( bamber.dims )
-
-#import matplotlib.pyplot as plt
-#plt.imshow(pri_data.filled(-1))
-#plt.show()
-
-print('Interpolating...')
-# loop through all data points
-for ii in range(0, bamber.ny):
-    for jj in range(0, bamber.nx):
-        # interpolate bamber point on massCon grid
-        indx = np.ravel_multi_index( (ii,jj),b_trans.dims )
-
-        nn_mc = b_trans.qi[indx]
-        nn_ii, nn_jj = np.unravel_index(nn_mc, massCon.dims)
-        
-        y_s = -1
-        x_s = -1
-        
-        # make sure inside priority grid (massCon)
-        if (b_trans.y_grid[ii,jj] < massCon.y_grid[0,0] or b_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
-                reason[ii,jj] = 1
-                continue
-            else:
-                tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid[ii,jj], bamber.y_grid[ii,jj], sec_data[ii,jj])
-                new_data[ii,jj] = td
-                reason[ii,jj] = -1
-                continue
-        
-        if (b_trans.x_grid[ii,jj] < massCon.x_grid[0,0] or b_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
-                reason[ii,jj] = 2
-                continue
-            else:
-                tx, ty, td = pyproj.transform(proj_eigen_gl04c, proj_epsg3413, bamber.x_grid[ii,jj], bamber.y_grid[ii,jj], sec_data[ii,jj])
-                new_data[ii,jj] = td
-                reason[ii,jj] = -1
-                continue
-
-        # find quadrent
-        if b_trans.y_grid[ii,jj] >= massCon.y_grid[nn_ii,nn_jj]:
-            y_s = +1
-        if b_trans.x_grid[ii,jj] >= massCon.x_grid[nn_ii,nn_jj]:
-            x_s = +1
-
-        # check for missing priority data!
-        missing_points, interp_dict = interp.check_missing(pri_data, (nn_ii, nn_jj), y_s, x_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, bamber.x_grid[ii,jj], bamber.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
-                reason[ii,jj] = -3
-            
-            else:
-                new_data[ii,jj] = td
-                reason[ii,jj] = -2
-                continue
-
-        else:
-            new_data.mask[ii,jj] = True
-            reason[ii,jj] = 3
-            continue
-
-        # interpolate!
-        alpha = ( b_trans.y_grid[ii,jj] - massCon.y_grid[nn_ii,nn_jj] )/(massCon.dy*y_s)
-        beta  = ( b_trans.x_grid[ii,jj] - massCon.x_grid[nn_ii,nn_jj] )/(massCon.dx*x_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+x_s) ]*w[1] \
-                         +interp_dict[ (nn_ii+y_s,nn_jj+x_s) ]*w[2] \
-                         +interp_dict[ (nn_ii+y_s,nn_jj    ) ]*w[3]
-
-import matplotlib.pyplot as plt
-plt.imshow(reason[::-1,:])
-plt.show()
-
-print('Writing data')
-test_var = nc_test.createVariable( 'interpData', 'f4', ('y','x',) )
-test_var[:,:] = new_data[:,:]
-
-test_reason = nc_test.createVariable('reason','i',('y','x',))
-test_reason[:,:] = reason[:,:]
-nc_test.close()
-print('And Done!')
diff --git a/tree_plot.py b/tree_plot.py
deleted file mode 100644
index 0955a18..0000000
--- a/tree_plot.py
+++ /dev/null
@@ -1,38 +0,0 @@
-import numpy as np
-from util.ncfunc import copy_atts, get_nc_file
-from util.projections import DataGrid
-
-
-nc_test = get_nc_file('test.nc','r')
-test = DataGrid()
-
-test.reason = nc_test.variables['reason']
-test.dims = test.reason[:,:].shape
-
-reason = np.ndarray(test.dims)
-reason[:,:] = test.reason[:,:]
-
-test.img = np.ones( (test.dims[0], test.dims[1], 3) )
-test.img *= 255
-
-# use secondary as interp data! 
-test.img[reason == -3, : ] = [255,0,0]
-# use secondary as new data! 
-test.img[reason == -1, : ] = [102,255,102]
-test.img[reason == -2, : ] = [102,255,102]
-# outside y range and no secondary data
-test.img[reason == 1, : ] = [255,102,255]
-# outside x range and no secondary data 
-test.img[reason == 2, : ] = [178,102,255]
-# inside but no secondary data
-test.img[reason == 3, : ] = [255,102,178]
-
-
-
-test.img /= 255
-
-import matplotlib.pyplot as plt
-plt.imshow(test.img[::-1,:,:])
-plt.show()
-
-nc_test.close