Add Figure.tilemap to plot XYZ tile maps (#2394)

Adding the tilemap method for plotting XYZ tile maps.
This is a wrapper around `pygmt.datasets.load_tile_map`
and `pygmt.Figure.grdimage`. Aliases from `grdimage` have
been copied over, and docstring for parameters from
`load_tile_map` have been copied over too.

* Add rioxarray to CI build matrix and include it as optional dependency

Let the Continuous Integration tests run with `rioxarray`,
include it in pyproject.toml and environment.yml, and document
it in `doc/install.rst` as an optional dependency.

* Reproject and plot map in lonlat coordinates by default

Given a region in lonlat coordinates, ensure that the output
figure is also plotted in lonlat instead of Web Mercator.

* Ensure that plot is clipped to bounding box region when no_clip is False

Pass the region to the grdimage call so that the plot extent is
the same as the bounding box region. Set `no_clip=True` to prevent
the clipping from happening (i.e., the plot will extend out from
the region of interest). Added a unit test checking results for
both no_clip True/False, and updated some previous baseline images
that have changed slightly.

* Web Mercator to Spherical Mercator and remove verbose statement

* Use rio.set_crs instead of rio.write_crs

So that the `spatial_ref` CF coordinate won't be set, which would
result in an int32 variable on Windows but int64 on Linux/macOS.

---------

Co-authored-by: Dongdong Tian <[email protected]>
Co-authored-by: Michael Grund <[email protected]>
Co-authored-by: Yvonne Fröhlich <[email protected]>

Author: 4 people

.github/workflows/ci_docs.yml

@@ -71,7 +71,7 @@ jobs:
       - name: Install dependencies
         run: |
           mamba install gmt=6.4.0 numpy pandas xarray netCDF4 packaging \
-                        build ipython make myst-parser contextily geopandas \
+                        build ipython make myst-parser contextily geopandas rioxarray \
                         sphinx sphinx-copybutton sphinx-design sphinx-gallery sphinx_rtd_theme
 
       # Show installed pkg information for postmortem diagnostic

.github/workflows/ci_tests.yaml

@@ -47,7 +47,7 @@ jobs:
             optional-packages: ''
           - python-version: '3.11'
             numpy-version: '1.24'
-            optional-packages: 'contextily geopandas ipython'
+            optional-packages: 'contextily geopandas ipython rioxarray'
     timeout-minutes: 30
     defaults:
       run:

.github/workflows/ci_tests_dev.yaml

@@ -101,8 +101,9 @@ jobs:
                         geopandas ghostscript libnetcdf hdf5 zlib curl pcre make
           pip install --pre --prefer-binary \
                         numpy pandas xarray netCDF4 packaging \
-                        build contextily dvc ipython 'pytest>=6.0' pytest-cov \
-                        pytest-doctestplus pytest-mpl sphinx-gallery
+                        build contextily dvc ipython rioxarray \
+                        'pytest>=6.0' pytest-cov pytest-doctestplus pytest-mpl \
+                        sphinx-gallery
 
       # Pull baseline image data from dvc remote (DAGsHub)
       - name: Pull baseline image data from dvc remote

.github/workflows/ci_tests_legacy.yaml

@@ -66,7 +66,7 @@ jobs:
         run: |
           mamba install gmt=${{ matrix.gmt_version }} numpy \
                         pandas xarray netCDF4 packaging \
-                        contextily geopandas ipython \
+                        contextily geopandas ipython rioxarray \
                         build dvc make 'pytest>=6.0' \
                         pytest-cov pytest-doctestplus pytest-mpl sphinx-gallery
 

ci/requirements/docs.yml

@@ -14,6 +14,7 @@ dependencies:
     # Optional dependencies
     - contextily
     - geopandas
+    - rioxarray
     # Development dependencies (general)
     - build
     - ipython

doc/api/index.rst

@@ -62,6 +62,7 @@ Plotting raster data
     Figure.grdimage
     Figure.grdview
     Figure.image
+    Figure.tilemap
 
 Configuring layout
 ~~~~~~~~~~~~~~~~~~

doc/conf.py

@@ -60,6 +60,7 @@
     "python": ("https://docs.python.org/3/", None),
     "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None),
     "rasterio": ("https://rasterio.readthedocs.io/en/stable/", None),
+    "rioxarray": ("https://corteva.github.io/rioxarray/stable/", None),
     "xarray": ("https://docs.xarray.dev/en/stable/", None),
     "xyzservices": ("https://xyzservices.readthedocs.io/en/stable", None),
 }

doc/install.rst

@@ -108,6 +108,7 @@ The following are optional dependencies:
 * `IPython <https://ipython.org>`__: For embedding the figures in Jupyter notebooks (recommended).
 * `Contextily <https://contextily.readthedocs.io>`__: For retrieving tile maps from the internet.
 * `GeoPandas <https://geopandas.org>`__: For using and plotting GeoDataFrame objects.
+* `RioXarray <https://corteva.github.io/rioxarray>`__: For saving multi-band rasters to GeoTIFFs.
 
 Installing GMT and other dependencies
 -------------------------------------

environment.yml

@@ -15,6 +15,7 @@ dependencies:
     - contextily
     - geopandas
     - ipython
+    - rioxarray
     # Development dependencies (general)
     - build
     - dvc

pygmt/datasets/tile_map.py

@@ -103,16 +103,16 @@ def load_tile_map(region, zoom="auto", source=None, lonlat=True, wait=0, max_ret
     Frozen({'band': 3, 'y': 256, 'x': 512})
     >>> raster.coords
     Coordinates:
-      * band     (band) uint8 0 1 2
-      * y        (y) float64 -7.081e-10 -7.858e+04 ... -1.996e+07 -2.004e+07
-      * x        (x) float64 -2.004e+07 -1.996e+07 ... 1.996e+07 2.004e+07
+      * band         (band) uint8 0 1 2
+      * y            (y) float64 -7.081e-10 -7.858e+04 ... -1.996e+07 ...
+      * x            (x) float64 -2.004e+07 -1.996e+07 ... 1.996e+07 2.004e+07
     """
     # pylint: disable=too-many-locals
     if contextily is None:
         raise ImportError(
             "Package `contextily` is required to be installed to use this function. "
             "Please use `pip install contextily` or "
-            "`conda install -c conda-forge contextily` "
+            "`mamba install -c conda-forge contextily` "
             "to install the package."
         )
 
@@ -147,6 +147,6 @@ def load_tile_map(region, zoom="auto", source=None, lonlat=True, wait=0, max_ret
 
     # If rioxarray is installed, set the coordinate reference system
     if hasattr(dataarray, "rio"):
-        dataarray = dataarray.rio.write_crs(input_crs="EPSG:3857")
+        dataarray = dataarray.rio.set_crs(input_crs="EPSG:3857")
 
     return dataarray

pygmt/figure.py

@@ -523,6 +523,7 @@ def _repr_html_(self):
         subplot,
         ternary,
         text,
+        tilemap,
         timestamp,
         velo,
         wiggle,

pygmt/src/__init__.py

@@ -51,6 +51,7 @@
 from pygmt.src.surface import surface
 from pygmt.src.ternary import ternary
 from pygmt.src.text import text_ as text  # "text" is an argument within "text_"
+from pygmt.src.tilemap import tilemap
 from pygmt.src.timestamp import timestamp
 from pygmt.src.triangulate import triangulate
 from pygmt.src.velo import velo

pygmt/src/tilemap.py

@@ -0,0 +1,155 @@
+"""
+tilemap - Plot XYZ tile maps.
+"""
+from pygmt.clib import Session
+from pygmt.datasets.tile_map import load_tile_map
+from pygmt.helpers import (
+    GMTTempFile,
+    build_arg_string,
+    fmt_docstring,
+    kwargs_to_strings,
+    use_alias,
+)
+
+try:
+    import rioxarray
+except ImportError:
+    rioxarray = None
+
+
+@fmt_docstring
+@use_alias(
+    B="frame",
+    E="dpi",
+    I="shading",
+    J="projection",
+    M="monochrome",
+    N="no_clip",
+    Q="nan_transparent",
+    # R="region",
+    V="verbose",
+    c="panel",
+    p="perspective",
+    t="transparency",
+)
+@kwargs_to_strings(c="sequence_comma", p="sequence")  # R="sequence",
+def tilemap(
+    self, region, zoom="auto", source=None, lonlat=True, wait=0, max_retries=2, **kwargs
+):
+    r"""
+    Plots an XYZ tile map.
+
+    This method loads XYZ tile maps from a tile server or local file using
+    :func:`pygmt.datasets.load_tile_map` into a georeferenced form, and plots
+    the tiles as a basemap or overlay using :meth:`pygmt.Figure.grdimage`.
+
+    **Note**: By default, standard web map tiles served in a Spherical Mercator
+    (EPSG:3857) Cartesian format will be reprojected to a geographic coordinate
+    reference system (OGC:WGS84) and plotted with longitude/latitude bounds
+    when ``lonlat=True``. If reprojection is not desired, please set
+    ``lonlat=False`` and provide Spherical Mercator (EPSG:3857) coordinates to
+    the ``region`` parameter.
+
+    {aliases}
+
+    Parameters
+    ----------
+    region : list
+        The bounding box of the map in the form of a list [*xmin*, *xmax*,
+        *ymin*, *ymax*]. These coordinates should be in longitude/latitude if
+        ``lonlat=True`` or Spherical Mercator (EPSG:3857) if ``lonlat=False``.
+
+    zoom : int or str
+        Optional. Level of detail. Higher levels (e.g. ``22``) mean a zoom
+        level closer to the Earth's surface, with more tiles covering a smaller
+        geographical area and thus more detail. Lower levels (e.g. ``0``) mean
+        a zoom level further from the Earth's surface, with less tiles covering
+        a larger geographical area and thus less detail [Default is
+        ``"auto"`` to automatically determine the zoom level based on the
+        bounding box region extent].
+
+        **Note**: The maximum possible zoom level may be smaller than ``22``,
+        and depends on what is supported by the chosen web tile provider
+        source.
+
+    source : xyzservices.TileProvider or str
+        Optional. The tile source: web tile provider or path to a local file.
+        Provide either:
+
+        - A web tile provider in the form of a
+          :class:`xyzservices.TileProvider` object. See
+          :doc:`Contextily providers <contextily:providers_deepdive>` for a
+          list of tile providers [Default is
+          ``xyzservices.providers.Stamen.Terrain``, i.e. Stamen Terrain web
+          tiles].
+        - A web tile provider in the form of a URL. The placeholders for the
+          XYZ in the URL need to be {{x}}, {{y}}, {{z}}, respectively. E.g.
+          ``https://{{s}}.tile.openstreetmap.org/{{z}}/{{x}}/{{y}}.png``.
+        - A local file path. The file is read with
+          :doc:`rasterio <rasterio:index>` and all bands are loaded into the
+          basemap. See
+          :doc:`contextily:working_with_local_files`.
+
+        IMPORTANT: Tiles are assumed to be in the Spherical Mercator projection
+        (EPSG:3857).
+
+    lonlat : bool
+        Optional. If ``False``, coordinates in ``region`` are assumed to be
+        Spherical Mercator as opposed to longitude/latitude [Default is
+        ``True``].
+
+    wait : int
+        Optional. If the tile API is rate-limited, the number of seconds to
+        wait between a failed request and the next try [Default is ``0``].
+
+    max_retries : int
+        Optional. Total number of rejected requests allowed before contextily
+        will stop trying to fetch more tiles from a rate-limited API [Default
+        is ``2``].
+
+    kwargs : dict
+        Extra keyword arguments to pass to :meth:`pygmt.Figure.grdimage`.
+
+    Raises
+    ------
+    ImportError
+        If ``rioxarray`` is not installed. Follow
+        :doc:`install instructions for rioxarray <rioxarray:installation>`,
+        (e.g. via ``pip install rioxarray``) before using this function.
+    """
+    kwargs = self._preprocess(**kwargs)  # pylint: disable=protected-access
+
+    if rioxarray is None:
+        raise ImportError(
+            "Package `rioxarray` is required to be installed to use this function. "
+            "Please use `pip install rioxarray` or "
+            "`mamba install -c conda-forge rioxarray` "
+            "to install the package."
+        )
+
+    raster = load_tile_map(
+        region=region,
+        zoom=zoom,
+        source=source,
+        lonlat=lonlat,
+        wait=wait,
+        max_retries=max_retries,
+    )
+
+    # Reproject raster from Spherical Mercator (EPSG:3857) to
+    # lonlat (OGC:CRS84) if bounding box region was provided in lonlat
+    if lonlat and raster.rio.crs == "EPSG:3857":
+        raster = raster.rio.reproject(dst_crs="OGC:CRS84")
+        raster.gmt.gtype = 1  # set to geographic type
+
+    # Only set region if no_clip is None or False, so that plot is clipped to
+    # exact bounding box region
+    if kwargs.get("N") in [None, False]:
+        kwargs["R"] = "/".join(str(coordinate) for coordinate in region)
+
+    with GMTTempFile(suffix=".tif") as tmpfile:
+        raster.rio.to_raster(raster_path=tmpfile.name)
+        with Session() as lib:
+            lib.call_module(
+                module="grdimage", args=build_arg_string(kwargs, infile=tmpfile.name)
+            )

pygmt/tests/baseline/test_tilemap_no_clip_False.png.dvc

@@ -0,0 +1,4 @@
+outs:
+- md5: 9317080021b0ce6f3b9ea6d17feece00
+  size: 23275
+  path: test_tilemap_no_clip_False.png

pygmt/tests/baseline/test_tilemap_no_clip_True.png.dvc

@@ -0,0 +1,4 @@
+outs:
+- md5: 83e6119b2351f9d472ca7e3cc45388c3
+  size: 60984
+  path: test_tilemap_no_clip_True.png

pygmt/tests/baseline/test_tilemap_ogc_wgs84.png.dvc

@@ -0,0 +1,4 @@
+outs:
+- md5: 3de0555d86aca49b92425c8d5272a934
+  size: 59286
+  path: test_tilemap_ogc_wgs84.png

pygmt/tests/baseline/test_tilemap_web_mercator.png.dvc

@@ -0,0 +1,4 @@
+outs:
+- md5: a76d9a9a1890d6b1345305eaea598bc3
+  size: 122195
+  path: test_tilemap_web_mercator.png

pygmt/tests/test_tilemap.py

@@ -0,0 +1,54 @@
+"""
+Tests Figure.tilemap.
+"""
+import pytest
+from pygmt import Figure
+
+contextily = pytest.importorskip("contextily")
+rioxarray = pytest.importorskip("rioxarray")
+
+
+@pytest.mark.mpl_image_compare
+def test_tilemap_web_mercator():
+    """
+    Create a tilemap plot in Spherical Mercator projection (EPSG:3857).
+    """
+    fig = Figure()
+    fig.tilemap(
+        region=[-20000000.0, 20000000.0, -20000000.0, 20000000.0],
+        zoom=0,
+        lonlat=False,
+        frame="afg",
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+def test_tilemap_ogc_wgs84():
+    """
+    Create a tilemap plot using longitude/latitude coordinates (OGC:WGS84),
+    centred on the international date line.
+    """
+    fig = Figure()
+    fig.tilemap(
+        region=[-180.0, 180.0, -90, 90], zoom=0, frame="afg", projection="R180/5c"
+    )
+    return fig
+
+
+@pytest.mark.mpl_image_compare
+@pytest.mark.parametrize("no_clip", [False, True])
+def test_tilemap_no_clip(no_clip):
+    """
+    Create a tilemap plot clipped to the Southern Hemisphere when no_clip is
+    False, but for the whole globe when no_clip is True.
+    """
+    fig = Figure()
+    fig.tilemap(
+        region=[-180.0, 180.0, -90, 0.6886],
+        zoom=0,
+        frame="afg",
+        projection="H180/5c",
+        no_clip=no_clip,
+    )
+    return fig

pyproject.toml

@@ -16,7 +16,7 @@ keywords = [
     "geophysics",
     "geospatial",
     "oceanography",
-    "seismology"
+    "seismology",
 ]
 classifiers = [
     "Development Status :: 4 - Beta",
@@ -36,15 +36,16 @@ dependencies = [
     "pandas",
     "xarray",
     "netCDF4",
-    "packaging"
+    "packaging",
 ]
 dynamic = ["version"]
 
 [project.optional-dependencies]
 all = [
     "contextily",
     "geopandas",
-    "ipython"
+    "ipython",
+    "rioxarray",
 ]
 
 [project.urls]