Add extent in polygon

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_in_polygon(crs=crs)
+        _i, _j = poly.exterior.xy
         return [np.min(_i), np.max(_i), np.min(_j), np.max(_j)]
 
+    def extent_in_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

@@ -723,6 +723,13 @@ def test_extent(self):
         assert_allclose([np.min(lon), np.min(lat)], [exgx[0], exgy[0]],
                         rtol=0.1)
 
+        p = g2.extent_in_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_in_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_in_polygon(), edgecolor='r',
+                    linewidth=2)
+    sm.visualize(ax=ax)
     return fig