Add reproject method

.travis.yml

@@ -27,9 +27,8 @@ install:
   - pip install codecov
   - pip install .
 
-
 script:
-  - coverage run --source pygridtools check_pygridtools.py --verbose --pep8 --mpl
+  - coverage run --source pygridtools check_pygridtools.py --verbose --strict
 
 after_success:
   - if [ ${COVERAGE} = true ]; then

pygridtools/core.py

@@ -214,6 +214,9 @@ class ModelGrid(object):
     ----------
     nodes_x, nodes_y : numpy.ndarray
         M-by-N arrays of node (vertex) coordinates for the grid.
+    crs : string
+        proj4-compliant coordinate reference system specification. See
+        http://geopandas.org/projections.html for more information.
 
     """
     def __init__(self, nodes_x, nodes_y, crs=None):
@@ -251,20 +254,12 @@ def nodes_y(self, value):
     @property
     def cells_x(self):
         """Array of cell centroid x-coordinates"""
-        xc = 0.25 * (
-            self.xn[1:, 1:] + self.xn[1:, :-1] +
-            self.xn[:-1, 1:] + self.xn[:-1, :-1]
-        )
-        return xc
+        return 0.25 * misc.padded_sum(self.xn)
 
     @property
     def cells_y(self):
         """Array of cell centroid y-coordinates"""
-        yc = 0.25 * (
-            self.yn[1:, 1:] + self.yn[1:, :-1] +
-            self.yn[:-1, 1:] + self.yn[:-1, :-1]
-        )
-        return yc
+        return 0.25 * misc.padded_sum(self.yn)
 
     @property
     def shape(self):
@@ -278,42 +273,42 @@ def cell_shape(self):
 
     @property
     def xn(self):
-        """Shortcut to x-coords of nodes"""
+        """ Shortcut to x-coords of nodes """
         return self.nodes_x
 
     @property
     def yn(self):
-        """Shortcut to y-coords of nodes"""
+        """ Shortcut to y-coords of nodes """
         return self.nodes_y
 
     @property
     def xc(self):
-        """ Shortcut to x-coords of cells/centroids"""
+        """ Shortcut to x-coords of cells/centroids """
         return self.cells_x
 
     @property
     def yc(self):
-        """ Shortcut to y-coords of cells/centroids"""
+        """ Shortcut to y-coords of cells/centroids """
         return self.cells_y
 
     @property
     def icells(self):
-        """ Number of rows of cells"""
+        """ Number of rows of cells """
         return self.cell_shape[1]
 
     @property
     def jcells(self):
-        """ Number of columns of cells"""
+        """ Number of columns of cells """
         return self.cell_shape[0]
 
     @property
     def inodes(self):
-        """Number of rows of nodes"""
+        """Number of rows of nodes """
         return self.shape[1]
 
     @property
     def jnodes(self):
-        """Number of columns of nodes"""
+        """Number of columns of nodes """
         return self.shape[0]
 
     @property
@@ -325,15 +320,17 @@ def cell_mask(self):
     def cell_mask(self, value):
         self._cell_mask = value
 
+    @property
+    def node_mask(self):
+        padded = numpy.pad(self.cell_mask, pad_width=1, mode='edge')
+        windowed = misc.padded_sum(padded.astype(int), window=1)
+        return (windowed == 4)
+
     @property
     def crs(self):
         """ Coordinate reference system for GIS data export """
         return self._crs
 
-    @crs.setter
-    def crs(self, value):
-        self._crs = value
-
     @property
     def domain(self):
         """ The optional domain used to generate the raw grid """
@@ -543,9 +540,10 @@ def extract(self, jstart=0, istart=0, jend=-1, iend=-1):
         return self.transform(extract, jstart=jstart, istart=istart, jend=jend, iend=iend)
 
     def copy(self):
+        """ Returns a deep copy of the current ModelGrid """
         return deepcopy(self)
 
-    def merge(self, other, how='vert', where='+', shift=0):
+    def merge(self, other, how='vert', where='+', shift=0, min_nodes=1):
         """
         Merge with another grid using pygridtools.misc.padded_stack.
 
@@ -581,6 +579,8 @@ def merge(self, other, how='vert', where='+', shift=0):
             axis other than the one being merged. In other words,
             vertically stacked arrays can be shifted horizontally,
             and horizontally stacked arrays can be shifted vertically.
+        min_nodes : int (default = 1)
+            Minimum number of masked nodes required to mask a cell.
 
         Returns
         -------
@@ -620,10 +620,13 @@ def merge(self, other, how='vert', where='+', shift=0):
 
         """
 
+        node_mask = merge(self.node_mask, other.node_mask, how=how, where=where, shift=shift)
+        cell_mask = misc.padded_sum(node_mask) >= min_nodes
+
         return ModelGrid(
             merge(self.nodes_x, other.nodes_x, how=how, where=where, shift=shift),
             merge(self.nodes_y, other.nodes_y, how=how, where=where, shift=shift)
-        ).update_cell_mask()
+        ).update_cell_mask(mask=cell_mask)
 
     def update_cell_mask(self, mask=None, merge_existing=True):
         """
@@ -694,16 +697,12 @@ def mask_nodes(self, polyverts, min_nodes=3, inside=False, use_existing=False,
 
         _node_mask = misc.mask_with_polygon(self.xn, self.yn, polyverts,
                                             inside=inside).astype(int)
-        cell_mask = (
-            _node_mask[1:, 1:] + _node_mask[:-1, :-1] +
-            _node_mask[:-1, 1:] + _node_mask[1:, :-1]
-        ) >= min_nodes
-        cell_mask = cell_mask.astype(bool)
+
+        cell_mask = (misc.padded_sum(_node_mask, window=1) >= min_nodes).astype(bool)
 
         return self.update_cell_mask(mask=cell_mask, merge_existing=use_existing)
 
     def mask_centroids(self, polyverts, inside=True, use_existing=True):
-
         """ Create mask for the cells of the ModelGrid with a polygon.
 
         Parameters
@@ -736,10 +735,9 @@ def mask_cells_with_polygon(self, polyverts, use_centroids=True, **kwargs):
         else:
             return self.mask_nodes(polyverts, **kwargs)
 
-    def plot_cells(self, engine='mpl', ax=None,
-                   usemask=True, cell_kws=None,
-                   domain_kws=None, extent_kws=None,
-                   showisland=True, island_kws=None):
+    def plot_cells(self, engine='mpl', ax=None, usemask=True, showisland=True,
+                   cell_kws=None, domain_kws=None, extent_kws=None,
+                   island_kws=None):
         """
         Creates a figure of the cells, boundary, domain, and islands.
 
@@ -773,13 +771,16 @@ def plot_cells(self, engine='mpl', ax=None,
                              mask=self.cell_mask, **cell_kws)
 
         if domain_kws is not None:
-            fig = viz.plot_domain(data=self.domain, engine=engine, ax=ax, **domain_kws)
+            fig = viz.plot_domain(data=self.domain, engine=engine,
+                                  ax=ax, **domain_kws)
 
         if extent_kws:
-            fig = viz.plot_boundaries(extent=self.extent, engine=engine, ax=ax, **extent_kws)
+            fig = viz.plot_boundaries(extent=self.extent, engine=engine,
+                                      ax=ax, **extent_kws)
 
         if island_kws:
-            fig = viz.plot_boundaries(islands=self.islands, engine=engine, ax=ax, **island_kws)
+            fig = viz.plot_boundaries(islands=self.islands, engine=engine,
+                                      ax=ax, **island_kws)
 
         return fig
 
@@ -838,6 +839,9 @@ def make_cols(top_level):
         northing = pandas.DataFrame(y, index=index, columns=northing_cols)
         return easting.join(northing)
 
+    def to_geodataframe(self, usemask=True, which='nodes'):
+        pass
+
     def to_coord_pairs(self, usemask=False, which='nodes'):
         """
         Converts a grid to a long array of coordinates pairs.
@@ -932,6 +936,13 @@ def to_gis(self, outputfile, usemask=True, which='cells',
     def to_gefdc(self, directory):
         return GEFDCWriter(self, directory)
 
+    def reproject(self, crs):
+        return (
+            self.to_geodataframe(usemask=True, which='nodes')
+                .to_crs(crs)
+                .pipe(ModelGrid.from_geodataframe)
+        )
+
     @classmethod
     def from_dataframe(cls, df, icol='ii', jcol='jj',
                        xcol='easting', ycol='northing'):
@@ -957,6 +968,10 @@ def from_dataframe(cls, df, icol='ii', jcol='jj',
         xtab = df.reset_index()[all_cols].set_index([icol, jcol]).unstack(level=icol)
         return cls(xtab[xcol], xtab[ycol]).update_cell_mask()
 
+    @classmethod
+    def from_geodataframe(cls, gdf, icol='ii', jcol='jj'):
+        pass
+
     @classmethod
     def from_gis(cls, gisfile, icol='ii', jcol='jj'):
         """

pygridtools/misc.py

@@ -273,6 +273,11 @@ def padded_stack(a, b, how='vert', where='+', shift=0, padval=numpy.nan):
     return stacked
 
 
+def padded_sum(padded, window=1):
+    return (padded[window:, window:] + padded[:-window, :-window] +
+            padded[:-window, window:] + padded[window:, :-window])
+
+
 def mask_with_polygon(x, y, polyverts, inside=True):
     """ Mask x-y arrays inside or outside a polygon
 

pygridtools/tests/test_core.py

@@ -464,14 +464,64 @@ def test_ModelGrid_split_ax0(mg, simple_nodes):
     nptest.assert_array_equal(mgbottom.nodes_y, yn[3:, :])
 
 
+def test_ModelGrid_node_mask(simple_nodes):
+    g = core.ModelGrid(*simple_nodes).update_cell_mask()
+    expected = numpy.array([
+        [0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 1, 1, 1, 1]
+    ]).astype(bool)
+    nptest.assert_array_equal(expected, g.node_mask)
+
+
 def test_ModelGrid_merge(g1, g2, simple_nodes):
     g3 = g1.merge(g2, how='horiz', where='+', shift=2)
-    g4 = core.ModelGrid(*simple_nodes)
+    g4 = core.ModelGrid(*simple_nodes).update_cell_mask()
 
     nptest.assert_array_equal(g3.xn, g4.xn)
     nptest.assert_array_equal(g3.xc, g4.xc)
 
 
+@pytest.mark.mpl_image_compare(baseline_dir=BASELINE_IMAGES, tolerance=15)
+def test_ModelGrid_merge_with_mask(simple_nodes):
+    mg1 = core.ModelGrid(*simple_nodes).update_cell_mask()
+    mg2 = (
+        mg1.transform_x(lambda x: x + 1)
+           .transform_y(lambda y: y + 5)
+           .update_cell_mask(mask=mg1.cell_mask)
+    )
+
+    merged = mg1.merge(mg2, where='+', shift=1, min_nodes=1)
+    expected = numpy.array([
+        [0, 0, 1, 1, 1, 1, 1],
+        [0, 0, 1, 1, 1, 1, 1],
+        [0, 0, 0, 0, 0, 0, 1],
+        [0, 0, 0, 0, 0, 0, 1],
+        [0, 0, 1, 1, 1, 1, 1],
+        [0, 0, 1, 1, 1, 1, 1],
+        [0, 0, 1, 1, 1, 1, 1],
+        [0, 0, 1, 1, 1, 1, 1],
+        [1, 0, 1, 1, 1, 1, 1],
+        [1, 0, 0, 1, 1, 1, 1],
+        [1, 0, 0, 1, 1, 1, 1],
+        [1, 0, 0, 0, 0, 0, 0],
+        [1, 0, 0, 0, 0, 0, 0],
+        [1, 0, 0, 1, 1, 1, 1],
+        [1, 0, 0, 1, 1, 1, 1],
+        [1, 0, 0, 1, 1, 1, 1],
+        [1, 0, 0, 1, 1, 1, 1]
+    ]).astype(bool)
+    nptest.assert_array_equal(merged.cell_mask, expected)
+    fig, artists = merged.plot_cells()
+    return fig
+
+
 def test_ModelGrid_insert_3_ax0(mg):
     known_xnodes = numpy.ma.masked_invalid(numpy.array([
         [1.0, 1.5, 2.0, nan, nan, nan, nan],

pygridtools/tests/test_misc.py

@@ -294,6 +294,35 @@ def test_padded_stack_errors(stackgrids, how, where):
                           how=how, where=where, shift=2)
 
 
+def test_padded_sum():
+    mask = numpy.array([
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+        [0, 0, 0, 0, 1, 1, 1, 1],
+    ])
+
+    result = misc.padded_sum(mask, window=1)
+    expected = numpy.array([
+        [0, 0, 0, 2, 4, 4, 4],
+        [0, 0, 0, 2, 4, 4, 4],
+        [0, 0, 0, 2, 4, 4, 4],
+        [0, 0, 0, 1, 2, 2, 2],
+        [0, 0, 0, 0, 0, 0, 0],
+        [0, 0, 0, 1, 2, 2, 2],
+        [0, 0, 0, 2, 4, 4, 4],
+        [0, 0, 0, 2, 4, 4, 4],
+        [0, 0, 0, 2, 4, 4, 4]
+    ])
+    nptest.assert_array_equal(result, expected)
+
+
 @pytest.mark.parametrize(('inside', 'expected'), [
     (True, numpy.array([
         [0, 0, 0, 0, 0], [0, 1, 1, 1, 0],