Add extent_as_polygon (#114)

Add extent_as_polygon

docs/whats-new.rst

@@ -11,6 +11,9 @@ v0.3.0 (unreleased)
 Enhancements
 ~~~~~~~~~~~~
 
+- new :py:func:`~Grid.extent_as_polygon` method, which creates a polygon
+  drawing the contours of a Grid.
+
 
 Bug fixes
 ~~~~~~~~~

salem/__init__.py

@@ -46,7 +46,7 @@ def _lazy_property(self):
 if not path.exists(download_dir):
     makedirs(download_dir)
 
-sample_data_gh_commit = 'a8de4369eb803ff1413b5fe7b5d6ac130bdbe2bb'
+sample_data_gh_commit = '63067a04762ecc22f105adc57aea1caa1a465ae1'
 sample_data_dir = path.join(cache_dir, 'salem-sample-data-' +
                             sample_data_gh_commit)
 

salem/gis.py

@@ -529,17 +529,37 @@ def extent_in_crs(self, crs=wgs84):
 
         # this is not so trivial
         # for optimisation we will transform the boundaries only
-        _i = np.hstack([np.arange(self.nx),
-                        np.ones(self.ny)*self.nx,
-                        np.arange(self.nx),
-                        np.zeros(self.ny)]).flatten()
-        _j = np.hstack([np.zeros(self.nx),
-                        np.arange(self.ny),
-                        np.ones(self.nx)*self.ny,
-                        np.arange(self.ny)]).flatten()
-        _i, _j = self.corner_grid.ij_to_crs(_i, _j, crs=crs)
+        poly = self.extent_as_polygon(crs=crs)
+        _i, _j = poly.exterior.xy
         return [np.min(_i), np.max(_i), np.min(_j), np.max(_j)]
 
+    def extent_as_polygon(self, crs=wgs84):
+        """Get the extent of the grid in a shapely.Polygon and desired crs.
+
+        Parameters
+        ----------
+        crs : crs
+            the target coordinate reference system.
+
+        Returns
+        -------
+        [left, right, bottom, top] boundaries of the grid.
+        """
+        from shapely.geometry import Polygon
+
+        # this is not so trivial
+        # for optimisation we will transform the boundaries only
+        _i = np.hstack([np.arange(self.nx+1),
+                        np.ones(self.ny+1)*self.nx,
+                        np.arange(self.nx+1)[::-1],
+                        np.zeros(self.ny+1)]).flatten()
+        _j = np.hstack([np.zeros(self.nx+1),
+                        np.arange(self.ny+1),
+                        np.ones(self.nx+1)*self.ny,
+                        np.arange(self.ny+1)[::-1]]).flatten()
+        _i, _j = self.corner_grid.ij_to_crs(_i, _j, crs=crs)
+        return Polygon(zip(_i, _j))
+
     def regrid(self, nx=None, ny=None, factor=1):
         """Make a copy of the grid with an updated spatial resolution.
 

salem/tests/test_gis.py

@@ -704,6 +704,7 @@ def test_map_gridded_data(self):
             self.assertTrue(odata.shape == (ny-1, nx-1))
             assert_allclose(exp_data, odata, atol=1e-03)
 
+    @requires_shapely
     def test_extent(self):
 
         # It should work exact same for any projection
@@ -723,6 +724,13 @@ def test_extent(self):
         assert_allclose([np.min(lon), np.min(lat)], [exgx[0], exgy[0]],
                         rtol=0.1)
 
+        p = g2.extent_as_polygon(crs=g2.proj)
+
+        assert p.is_valid
+        x, y = p.exterior.coords.xy
+        assert_allclose([np.min(x), np.max(x), np.min(y), np.max(y)],
+                        g2.extent)
+
     def test_simple_dataset(self):
         # see if with is working
         with SimpleNcDataSet(get_demo_file('dem_wgs84.nc')) as nc:

salem/tests/test_graphics.py

@@ -277,6 +277,13 @@ def test_contourf():
 
     c.set_lonlat_contours(interval=0.5)
 
+    # Add a geometry for fun
+    gs = g.to_dict()
+    gs['nxny'] = (1, 2)
+    gs['x0y0'] = (0, 2)
+    gs = Grid.from_dict(gs)
+    c.set_geometry(gs.extent_as_polygon(), edgecolor='r', linewidth=2)
+
     fig, ax = plt.subplots(1)
     c.visualize(ax=ax)
     fig.tight_layout()
@@ -287,7 +294,6 @@ def test_contourf():
 
     return fig
 
-
 @requires_matplotlib
 @pytest.mark.mpl_image_compare(baseline_dir=baseline_dir)
 def test_merca_map():
@@ -657,7 +663,11 @@ def test_lookup_transform():
     dse = open_xr_dataset(get_demo_file('era_interim_tibet.nc'))
     out = dse.salem.lookup_transform(dsw.T2C.isel(time=0), method=len)
     fig, ax = plt.subplots(1, 1)
-    out.salem.quick_map(ax=ax)
+    sm = out.salem.get_map()
+    sm.set_data(out)
+    sm.set_geometry(dsw.salem.grid.extent_as_polygon(), edgecolor='r',
+                    linewidth=2)
+    sm.visualize(ax=ax)
     return fig