new transform method added to Map (#33)

docs/api.rst

@@ -94,6 +94,9 @@ Map methods
     Map.set_shapefile
     Map.set_text
     Map.set_topography
+    Map.transform
+    Map.visualize
+    Map.plot
 
 Input/output
 ============

docs/examples/plot_hydrosheds.py

@@ -0,0 +1,58 @@
+# -*- coding: utf-8 -*-
+"""
+===================
+Plotting shapefiles
+===================
+
+Put some colors on shapefiles
+
+In this script, we use data from the `HydroSHEDS <http://www.hydrosheds.org/>`_
+database to illustrate some functionalities of salem Maps.
+"""
+
+import salem
+import matplotlib.pyplot as plt
+
+shpf = salem.get_demo_file('Lev_09_MAIN_BAS_4099000881.shp')
+gdf = salem.read_shapefile(shpf)
+
+# Get the google map which encompasses all geometries
+g = salem.GoogleVisibleMap(x=[gdf.min_x.min(), gdf.max_x.max()],
+                           y=[gdf.min_y.min(), gdf.max_y.max()],
+                           maptype='satellite', size_x=400, size_y=400)
+ggl_img = g.get_vardata()
+
+# Get the polygon of the last sink (i.e. the lake)
+gds_0 = gdf.loc[gdf.HYBAS_ID == gdf.iloc[0].MAIN_BAS]
+
+# Get each level draining into the lake, then into the last level, and so on
+gds = []
+prev_id = [gdf.iloc[0].MAIN_BAS]
+while True:
+    gd = gdf.loc[gdf.NEXT_DOWN.isin(prev_id)]
+    if len(gd) == 0:
+        break
+    gds.append(gd)
+    prev_id = gd.HYBAS_ID.unique()
+
+# make a map of the same size as the image
+sm = salem.Map(g.grid, factor=1)
+sm.set_rgb(ggl_img)  # add the background rgb image
+sm.set_shapefile(gds_0, linewidth=3)  # add the lake outline
+# and add all the draining basins
+cmap = plt.get_cmap('Blues')
+for i, gd in enumerate(gds):
+    # here we use a trick. set_shapefile uses PatchCollections internally,
+    # which is fast but does not support legend labels.
+    # so we use set_geometry instead:
+    for g, geo in enumerate(gd.geometry):
+        # we don't want more than one label per level
+        label = 'Level {:02d}'.format(i+1) if g == 0 else None
+        sm.set_geometry(geo, facecolor=cmap(i/(len(gds)-1)),
+                        alpha=0.8, label=label)
+# plot!
+f, ax = plt.subplots(figsize=(6, 4))
+ax.set_position([0.05, 0.06, 0.7, 0.9])
+sm.visualize(ax=ax)
+plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+plt.show()

docs/whats-new.rst

@@ -4,6 +4,22 @@ What's New
 ==========
 
 
+v0.2.X (Unreleased)
+-------------------
+
+Enhancements
+~~~~~~~~~~~~
+
+- new ``Map.transform()`` method to make overplotting easier (still
+  a draft)
+
+Bug fixes
+~~~~~~~~~
+
+- ``grid.transform()`` now works with non numpy array type
+- ``transform_geopandas()`` won't do inplace operation per default anymore
+
+
 v0.2.0 (08 November 2016)
 -------------------------
 

salem/datasets.py

@@ -223,7 +223,8 @@ def set_roi(self, shape=None, geometry=None, crs=wgs84, grid=None,
         # Several cases
         if shape is not None:
             gdf = sio.read_shapefile(shape)
-            gis.transform_geopandas(gdf, to_crs=ogrid.corner_grid)
+            gis.transform_geopandas(gdf, to_crs=ogrid.corner_grid,
+                                    inplace=True)
             with rasterio.Env():
                 mask = features.rasterize(gdf.geometry, out=mask)
         if geometry is not None:

salem/gis.py

@@ -540,6 +540,8 @@ def transform(self, x, y, z=None, crs=wgs84, nearest=False, maskout=False):
         (i, j) coordinates of the points in the local grid.
         """
 
+        x, y = np.ma.array(x), np.ma.array(y)
+
         # First to local proj
         _crs = check_crs(crs)
         if isinstance(_crs, pyproj.Proj):
@@ -550,8 +552,8 @@ def transform(self, x, y, z=None, crs=wgs84, nearest=False, maskout=False):
             raise ValueError('crs not understood')
 
         # Then to local grid
-        x = (np.ma.array(x) - self.x0) / self.dx
-        y = (np.ma.array(y) - self.y0) / self.dy
+        x = (x - self.x0) / self.dx
+        y = (y - self.y0) / self.dy
 
         # See if we need to round
         if nearest:
@@ -735,11 +737,14 @@ def region_of_interest(self, shape=None, geometry=None, grid=None,
         # Several cases
         if shape is not None:
             import pandas as pd
+            inplace = False
             if not isinstance(shape, pd.DataFrame):
                 from salem.sio import read_shapefile
                 shape = read_shapefile(shape)
+                inplace = True
             # corner grid is needed for rasterio
-            transform_geopandas(shape, to_crs=self.corner_grid)
+            shape = transform_geopandas(shape, to_crs=self.corner_grid,
+                                        inplace=inplace)
             import rasterio
             with rasterio.Env():
                 mask = rasterio.features.rasterize(shape.geometry, out=mask)
@@ -853,7 +858,7 @@ def transform_geometry(geom, crs=wgs84, to_crs=wgs84):
     return shapely_transform(project, geom)
 
 
-def transform_geopandas(gdf, to_crs=wgs84, inplace=True):
+def transform_geopandas(gdf, to_crs=wgs84, inplace=False):
     """Reprojects a geopandas dataframe.
 
     Parameters
@@ -863,7 +868,7 @@ def transform_geopandas(gdf, to_crs=wgs84, inplace=True):
     to_crs : crs
         the crs into which the dataframe must be transformed
     inplace : bool
-        the original dataframe will be overwritten (default)
+        the original dataframe will be overwritten (default: False)
 
     Returns
     -------

salem/graphics.py

@@ -11,8 +11,6 @@
 import copy
 # External libs
 import numpy as np
-from numpy import ma
-
 
 try:
     from scipy.misc import imresize
@@ -32,6 +30,8 @@
     from matplotlib.collections import PatchCollection, LineCollection
     from shapely.geometry import MultiPoint
     from descartes.patch import PolygonPatch
+    from matplotlib.transforms import Transform as MPLTranform
+    import matplotlib.path as mpath
 except ImportError:
     class d1():
         def __init__(self):
@@ -40,6 +40,7 @@ class d2():
             self.colors = d2
             self.colors.BoundaryNorm = object
     mpl = d1()
+    MPLTranform = object
 
 from salem import utils, gis, sio, Grid, wgs84, cache_dir, GeoTiff
 
@@ -90,7 +91,7 @@ def __init__(self, boundaries, ncolors, extend='neither'):
 
     def __call__(self, value):
         xx, is_scalar = self.process_value(value)
-        mask = ma.getmaskarray(xx)
+        mask = np.ma.getmaskarray(xx)
         xx = np.atleast_1d(xx.filled(self.vmax + 1))
         iret = np.zeros(xx.shape, dtype=np.int16)
         for i, b in enumerate(self._b):
@@ -100,7 +101,7 @@ def __call__(self, value):
             iret = (iret * scalefac).astype(np.int16)
         iret[xx < self.vmin] = -1
         iret[xx >= self.vmax] = self.Ncmap
-        ret = ma.array(iret, mask=mask)
+        ret = np.ma.array(iret, mask=mask)
         if is_scalar:
             ret = int(ret[0])  # assume python scalar
         return ret
@@ -642,7 +643,7 @@ def set_shapefile(self, shape=None, countries=False, oceans=False,
             kwargs.setdefault('colors', (0.08984375, 0.65625, 0.8515625))
             return self.set_shapefile(shapefiles['rivers'], **kwargs)
         if countries:
-            kwargs.setdefault('zorder', 98)
+            kwargs.setdefault('zorder', 50)
             return self.set_shapefile(shapefiles['world_borders'], **kwargs)
 
         # Reset?
@@ -784,7 +785,7 @@ def set_lonlat_contours(self, interval=None, xinterval=None,
         # Done
         kwargs.setdefault('colors', 'gray')
         kwargs.setdefault('linestyles', 'dashed')
-        kwargs.setdefault('zorder', 99)
+        kwargs.setdefault('zorder', 51)
         self.ll_contour_kw = kwargs
 
     def _shading_base(self, slope=None, relief_factor=0.7):
@@ -892,6 +893,30 @@ def set_rgb(self, img=None, crs=None, interp='nearest',
             out.append(self._check_data(img[..., i], crs=crs, interp=interp))
         self._rgb = np.dstack(out)
 
+    def transform(self, crs=wgs84, ax=None):
+        """Get a matplotlib transform object for a given reference system
+
+        Parameters
+        ----------
+        crs : coordinate reference system
+            a Grid or a Proj, basically. If a grid is given, the grid's proj
+            will be used.
+
+        Returns
+        -------
+        a matplotlib.transforms.Transform instance
+        """
+        try:
+            crs = crs.salem.grid
+        except:
+            pass
+        try:
+            crs = crs.proj
+        except:
+            pass
+        return _SalemTransform(target_grid=self.grid,
+                               source_crs=crs, ax=ax)
+
     def to_rgb(self):
         """Transform the data to a RGB image and add topographical shading."""
 
@@ -967,6 +992,7 @@ def plot(self, ax):
                 kwargs.setdefault('color', 'k')
                 ax.plot(a[:, 0], a[:, 1], **kwargs)
             if g.type == 'Point':
+                kwargs.setdefault('zorder', 99)
                 kwargs.setdefault('marker', 'o')
                 kwargs.setdefault('s', 60)
                 kwargs.setdefault('facecolor', 'w')
@@ -1014,3 +1040,40 @@ def plot_polygon(ax, poly, edgecolor='black', **kwargs):
     for p in poly.interiors:
         x, y = zip(*p.coords)
         ax.plot(x, y, color=edgecolor)
+
+
+class _SalemTransform(MPLTranform):
+    """
+    A transform class for mpl axes using Grids.
+    """
+
+    input_dims = 2
+    output_dims = 2
+    is_separable = False
+    has_inverse = False
+
+    def __init__(self, target_grid=None, source_crs=None, ax=None):
+        """ Instanciate.
+
+        Parameters
+        ----------
+        target_grid : salem.Grid
+            typically, the map grid
+        source_grid
+        """
+        self.source_crs = source_crs
+        self.target_grid = target_grid
+        self.ax = ax
+        MPLTranform.__init__(self)
+
+    def transform_non_affine(self, xy):
+        xx, yy = xy[:, 0:1], xy[:, 1:2]
+        xx, yy = self.target_grid.transform(xx, yy, crs=self.source_crs)
+        ax = plt.gca() if self.ax is None else self.ax
+        return ax.transData.transform(np.concatenate((xx, yy), 1))
+
+    def transform_path_non_affine(self, path):
+        raise NotImplementedError('path transforms not working yet')
+
+    def inverted(self):
+        raise NotImplementedError('inverted not working yet')

salem/tests/test_gis.py

@@ -814,12 +814,12 @@ def test_shape(self):
 
         so = read_shapefile(get_demo_file('Hintereisferner.shp'))
         sref = read_shapefile(get_demo_file('Hintereisferner_UTM.shp'))
-        st = gis.transform_geopandas(so, to_crs=sref.crs, inplace=False)
+        st = gis.transform_geopandas(so, to_crs=sref.crs)
         self.assertFalse(st is so)
         assert_allclose(st.geometry[0].exterior.coords,
                                    sref.geometry[0].exterior.coords)
 
-        sti = gis.transform_geopandas(so, to_crs=sref.crs)
+        sti = gis.transform_geopandas(so, to_crs=sref.crs, inplace=True)
         self.assertTrue(sti is so)
         assert_allclose(so.geometry[0].exterior.coords,
                                    sref.geometry[0].exterior.coords)
@@ -828,7 +828,7 @@ def test_shape(self):
 
         g = Grid(nxny=(1, 1), dxdy=(1, 1), ll_corner=(10., 46.), proj=wgs84)
         so = read_shapefile(get_demo_file('Hintereisferner.shp'))
-        st = gis.transform_geopandas(so, to_crs=g, inplace=False)
+        st = gis.transform_geopandas(so, to_crs=g)
 
         ref = np.array(so.geometry[0].exterior.coords)
         ref = ref - np.floor(ref)

salem/tests/test_graphics.py

@@ -654,8 +654,7 @@ def test_cartopy():
 
     ax = plt.subplot(3, 2, 3)
     smap = ds.salem.quick_map(ax=ax)
-    x, y = smap.grid.transform(ds.lon.values, ds.lat.values)
-    ax.scatter(x, y)
+    ax.scatter(ds.lon, ds.lat, transform=smap.transform(ax=ax))
 
     p = ds.salem.cartopy()
     ax = plt.subplot(3, 2, 4, projection=p)
@@ -668,8 +667,7 @@ def test_cartopy():
 
     ax = plt.subplot(3, 2, 5)
     smap = ds.salem.quick_map(ax=ax, factor=1)
-    x, y = smap.grid.transform(ds.lon.values, ds.lat.values)
-    ax.scatter(x, y)
+    ax.scatter(ds.lon, ds.lat, transform=smap.transform(ax=ax))
 
     p = ds.salem.cartopy()
     ax = plt.subplot(3, 2, 6, projection=p)