new example

README.rst

@@ -4,14 +4,16 @@
 Salem
 =====
 
+.. image:: https://badge.fury.io/py/salem.svg
+    :target: https://badge.fury.io/py/salem
 .. image:: https://travis-ci.org/fmaussion/salem.svg?branch=master
     :target: https://travis-ci.org/fmaussion/salem
 .. image:: https://coveralls.io/repos/fmaussion/salem/badge.svg?branch=master&service=github
-  :target: https://coveralls.io/github/fmaussion/salem?branch=master
+    :target: https://coveralls.io/github/fmaussion/salem?branch=master
 .. image:: https://readthedocs.org/projects/salem/badge/?version=stable
-  :target: http://salem.readthedocs.io/en/stable/?badge=stable
+    :target: http://salem.readthedocs.io/en/stable/?badge=stable
 .. image:: https://zenodo.org/badge/42607422.svg
-   :target: https://zenodo.org/badge/latestdoi/42607422
+    :target: https://zenodo.org/badge/latestdoi/42607422
 
 
 Salem is a `cat`_. Salem is also a small library to do some geoscientific data

docs/examples/README.txt

@@ -3,7 +3,7 @@
 Recipes
 =======
 
-Plot examples
--------------
+Plotting
+--------
 
-General-purpose examples
+How to make nice looking plots.

docs/examples/plot_googlestatic.py

@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+"""
+===========================
+Plot on a google static map
+===========================
+
+In this script, we use `motionless <https://github.com/ryancox/motionless>`_
+to download an image from the `google static map API <http://code.google.com/apis/maps/documentation/staticmaps/>`_
+and plot it on a :py:class:`~salem.Map`. We then add information to the map
+such as a glacier outline (from the `RGI <http://www.glims.org/RGI/>`_),
+and ground penetrating radar measurements (GPR, from the
+`GlaThiDa <http://www.gtn-g.ch/data_catalogue_glathida/>`_ database).
+
+
+The GPR measurements were realized in 1997, the glacier outlines are from
+2003, and the map background is from 2016. This illustrates the retreat of
+the Kesselwandferner glacier.
+
+"""
+
+import numpy as np
+import pandas as pd
+import salem
+from salem import get_demo_file, DataLevels, GoogleVisibleMap, Map
+import matplotlib.pyplot as plt
+
+# prepare the figure
+f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
+
+# read the shapefile and use its extent to define a ideally sized map
+shp = salem.read_shapefile(get_demo_file('rgi_kesselwand.shp'))
+g = GoogleVisibleMap(x=[shp.min_x, shp.max_x], y=[shp.min_y, shp.max_y],
+                     maptype='satellite')  # try out: 'terrain'
+
+# the google static image is a standard rgb image
+ggl_img = g.get_vardata()
+ax1.imshow(ggl_img)
+ax1.set_title('Google static map')
+
+# make a map of the same size as the image
+sm = Map(g.grid, factor=1)
+sm.set_shapefile(shp)  # add the glacier outlines
+sm.set_rgb(ggl_img)  # add the background rgb image
+sm.visualize(ax=ax2)  # plot it
+ax2.set_title('GPR measurements')
+
+# read the point GPR data and add them to the plot
+df = pd.read_csv(get_demo_file('gtd_ttt_kesselwand.csv'))
+dl = DataLevels(df.THICKNESS, levels=np.arange(10, 201, 10), extend='both')
+x, y = sm.grid.transform(df.POINT_LON.values, df.POINT_LAT.values)
+ax2.scatter(x, y, color=dl.to_rgb(), s=50, edgecolors='k', linewidths=1)
+dl.append_colorbar(ax2, label='Ice thickness (m)')
+
+# make it nice
+plt.tight_layout()
+plt.show()

docs/examples/plot_topo_shading.py

@@ -1,24 +1,44 @@
 # -*- coding: utf-8 -*-
 """
-==============
-Add topography
-==============
+===================
+Topographic shading
+===================
 
-Here we add topography to the plot.
+How to add topographic shading to a plot.
 
 """
 
 from salem import mercator_grid, Map, get_demo_file
 import matplotlib.pyplot as plt
 
-grid = mercator_grid(center_ll=(10.76, 46.798444),
-                     extent=(10000, 7000))
-c = Map(grid, countries=False)
-c.set_lonlat_contours(interval=10)
-c.set_shapefile(get_demo_file('Hintereisferner_UTM.shp'))
-c.set_topography(get_demo_file('hef_srtm.tif'))
+# prepare the figure
 
-fig, ax = plt.subplots(1, 1)
-c.visualize(ax=ax, addcbar=False, title='Default: spline deg 3')
+f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(8, 7))
+
+# map extent
+grid = mercator_grid(center_ll=(10.76, 46.79), extent=(18000, 14000))
+sm = Map(grid, countries=False)
+sm.set_lonlat_contours(interval=0)
+
+# add topography
+fpath = get_demo_file('hef_srtm.tif')
+sm.set_topography(fpath)
+sm.visualize(ax=ax1, addcbar=False, title='relief_factor=0.7 (default)')
+
+# stronger shading
+sm.set_topography(fpath, relief_factor=1.4)
+sm.visualize(ax=ax2, addcbar=False, title='relief_factor=1.4')
+
+# add color shading
+z = sm.set_topography(fpath)
+sm.set_data(z)
+sm.set_cmap('topo')
+sm.visualize(ax=ax3, title='Color with shading', addcbar=False)
+
+# color without topo shading
+sm.set_topography()
+sm.visualize(ax=ax4, title='Color without shading', addcbar=False)
+
+# make it nice
 plt.tight_layout()
 plt.show()

docs/whats-new.rst

@@ -8,12 +8,14 @@ v0.1.1 (Unreleased)
 Enhancements
 ~~~~~~~~~~~~
 
-- General documentation enhancements
+- General doc improvements
 - New ``ds`` keyword to the accessors ``subset()`` and ``roi()`` methods
 
 Bug fixes
 ~~~~~~~~~
 
+- Natural Earth file `lr` (low-res) now shipped with sample-data, `mr` and `hr`
+  can still be downloaded if needed
 - Remove use to deprecated ``rasterio.drivers()``
 - GeoNetCDF files without time variable should now open without error
 

salem/datasets.py

@@ -470,14 +470,32 @@ def __init__(self, file, grid=None, time=None):
 
 
 class GoogleCenterMap(GeoDataset):
-    """Google Static Maps (needs motionless)."""
+    """Google static map centered on a point.
+
+    Needs motionless.
+    """
 
     def __init__(self, center_ll=(11.38, 47.26), size_x=640, size_y=640,
                  zoom=12, maptype='satellite', use_cache=True, **kwargs):
-        """Open the file.
+        """Initialize
 
         Parameters
         ----------
+        center_ll : tuple
+          tuple of lon, lat center of the map
+        size_x : int
+          image size
+        size_y : int
+          image size
+        zoom int
+          google zoom level (https://developers.google.com/maps/documentation/
+          static-maps/intro#Zoomlevels)
+        maptype : str, default: 'satellite'
+          'roadmap', 'satellite', 'hybrid', 'terrain'
+        use_cache : bool, default: True
+          store the downloaded image in the cache to avoid future downloads
+        kwargs : **
+          any keyword accepted by motionless.CenterMap (e.g. `key` for the API)
         """
 
         if 'scale' in kwargs:
@@ -514,19 +532,44 @@ def get_vardata(self, var_id=0):
 
 
 class GoogleVisibleMap(GoogleCenterMap):
-    """Google Static Maps (needs motionless)."""
+    """Google static map automatically sized and zoomed to a selected region.
+
+    It's usually more practical to use than GoogleCenterMap.
+    """
 
-    def __init__(self, x, y, src=wgs84, size_x=640, size_y=640, **kwargs):
+    def __init__(self, x, y, crs=wgs84, size_x=640, size_y=640,
+                 maptype='satellite', use_cache=True, **kwargs):
+        """Initialize
+
+        Parameters
+        ----------
+        x : array
+          x coordinates of the points to include on the map
+        y : array
+          y coordinates of the points to include on the map
+        crs : proj or Grid
+          coordinate reference system of x, y
+        size_x : int
+          image size
+        size_y : int
+          image size
+        maptype : str, default: 'satellite'
+          'roadmap', 'satellite', 'hybrid', 'terrain'
+        use_cache : bool, default: True
+          store the downloaded image in the cache to avoid future downloads
+        kwargs : **
+          any keyword accepted by motionless.CenterMap (e.g. `key` for the API)
+        """
 
-        if 'zoom' in kwargs:
-            raise ValueError('zoom kwarg not accepted.')
+        if 'zoom' in kwargs or 'center_ll' in kwargs:
+            raise ValueError('incompatible kwargs.')
 
         # Transform to lonlat
-        src = gis.check_crs(src)
-        if isinstance(src, pyproj.Proj):
-            lon, lat = gis.transform_proj(src, wgs84, x, y)
-        elif isinstance(src, Grid):
-            lon, lat = src.ij_to_crs(x, y, crs=wgs84)
+        crs = gis.check_crs(crs)
+        if isinstance(crs, pyproj.Proj):
+            lon, lat = gis.transform_proj(crs, wgs84, x, y)
+        elif isinstance(crs, Grid):
+            lon, lat = crs.ij_to_crs(x, y, crs=wgs84)
         else:
             raise NotImplementedError()
 
@@ -535,12 +578,13 @@ def __init__(self, x, y, src=wgs84, size_x=640, size_y=640, **kwargs):
         zoom = 20
         while zoom >= 0:
             grid = gis.googlestatic_mercator_grid(center_ll=mc, nx=size_x,
-                                                    ny=size_y, zoom=zoom)
+                                                  ny=size_y, zoom=zoom)
             dx, dy = grid.transform(lon, lat, maskout=True)
             if np.any(dx.mask):
                 zoom -= 1
             else:
                 break
 
         GoogleCenterMap.__init__(self, center_ll=mc, size_x=size_x,
-                                 size_y=size_y, zoom=zoom, **kwargs)
+                                 size_y=size_y, zoom=zoom, maptype=maptype,
+                                 use_cache=use_cache, **kwargs)

salem/graphics.py

@@ -387,20 +387,32 @@ class Map(DataLevels):
     regional maps.
     """
 
-    def __init__(self, grid, nx=500, ny=None, countries=True, **kwargs):
+    def __init__(self, grid, nx=500, ny=None, factor=None,
+                 countries=True, **kwargs):
         """Make a new map.
 
         Parameters
         ----------
-        grid: a salem.Grid instance defining the map
-        nx: x resolution (in pixels) of the map
-        ny: y resolution (in pixels) of the map (ignored if nx is set)
-        countries: automatically add country borders to the map (you can do
-        it later with a call to set_shapefile)
-        kwards: all keywords accepted by DataLevels
+        grid : salem.Grid
+          the grid defining the map
+        nx : int
+          x resolution (in pixels) of the map
+        ny : int
+          y resolution (in pixels) of the map (ignored if nx is set)
+        factor : float
+          shortcut to keep the same image resolution as the grid
+        countries : bool
+          automatically add country borders to the map (you can do
+          it later with a call to set_shapefile)
+        kwargs: **
+          all keywords accepted by DataLevels
         """
 
-        self.grid = grid.center_grid.regrid(nx=nx, ny=ny)
+        if factor is not None:
+            nx = None
+            ny = None
+
+        self.grid = grid.center_grid.regrid(nx=nx, ny=ny, factor=factor)
         self.origin = 'lower' if self.grid.order == 'll' else 'upper'
 
         DataLevels.__init__(self, **kwargs)
@@ -698,9 +710,11 @@ def set_lonlat_contours(self, interval=None, xinterval=None,
                             **kwargs):
         """Add longitude and latitude contours to the map.
 
+        Calling it with interval=0. remove all conours.
+
         Parameters
         ----------
-        interval: interval (in degrees) between the contours (same for lon
+        interval : interval (in degrees) between the contours (same for lon
         and lat)
         xinterval: set a different interval for lons
         yinterval: set a different interval for lats
@@ -716,14 +730,22 @@ def set_lonlat_contours(self, interval=None, xinterval=None,
         if yinterval is None:
             yinterval = interval
 
-        # Change XY into interval coordinates, and back after rounding
-        xx, yy = self.grid.pixcorner_ll_coordinates
-        _xx = xx / xinterval
-        _yy = yy / yinterval
-        mm_x = [np.ceil(np.min(_xx)), np.floor(np.max(_xx))]
-        mm_y = [np.ceil(np.min(_yy)), np.floor(np.max(_yy))]
-        self.xtick_levs = (mm_x[0] + np.arange(mm_x[1]-mm_x[0]+1)) * xinterval
-        self.ytick_levs = (mm_y[0] + np.arange(mm_y[1]-mm_y[0]+1)) * yinterval
+        if xinterval == 0:
+            # no contour
+            self.xtick_levs = []
+            self.ytick_levs = []
+            add_tick_labels = False
+        else:
+            # Change XY into interval coordinates, and back after rounding
+            xx, yy = self.grid.pixcorner_ll_coordinates
+            _xx = xx / xinterval
+            _yy = yy / yinterval
+            mm_x = [np.ceil(np.min(_xx)), np.floor(np.max(_xx))]
+            mm_y = [np.ceil(np.min(_yy)), np.floor(np.max(_yy))]
+            self.xtick_levs = (mm_x[0] + np.arange(mm_x[1]-mm_x[0]+1)) * \
+                              xinterval
+            self.ytick_levs = (mm_y[0] + np.arange(mm_y[1]-mm_y[0]+1)) * \
+                              yinterval
 
         # Decide on float format
         d = np.array(['4', '3', '2', '1', '0'])
@@ -799,6 +821,8 @@ def set_topography(self, topo=None, crs=None, relief_factor=0.7, **kwargs):
 
         if topo is None:
             self._shading_base()
+            return
+
         kwargs.setdefault('interp', 'spline')
 
         if isinstance(topo, six.string_types):
@@ -844,7 +868,8 @@ def set_rgb(self, img=None, crs=None, interp='nearest',
         interp : str, default 'nearest'
             'nearest', 'linear', or 'spline'
         natural_earth : str
-           'lr' or 'hr' (low res or high res) natural earth background img
+           'lr', 'mr' or 'hr' (low res, medium or high res)
+           natural earth background img
         """
 
         if natural_earth is not None:

salem/tests/test_datasets.py

@@ -468,7 +468,7 @@ def test_visible(self):
         fw = get_demo_file('wrf_tip_d1.nc')
         d = GeoNetcdf(fw)
         i, j = d.grid.ij_coordinates
-        g = GoogleVisibleMap(x=i, y=j, src=d.grid, size_x=500, size_y=500)
+        g = GoogleVisibleMap(x=i, y=j, crs=d.grid, size_x=500, size_y=500)
         img = g.get_vardata()
         mask = g.grid.map_gridded_data(i*0+1, d.grid)