Add function to load raster tile maps using contextily

.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 geopandas \
+                        build ipython make myst-parser contextily geopandas \
                         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: 'geopandas ipython'
+            optional-packages: 'contextily geopandas ipython'
     timeout-minutes: 30
     defaults:
       run:

.github/workflows/ci_tests_dev.yaml

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

.github/workflows/ci_tests_legacy.yaml

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

ci/requirements/docs.yml

@@ -12,6 +12,7 @@ dependencies:
     - netCDF4
     - packaging
     # Optional dependencies
+    - contextily
     - geopandas
     # Development dependencies (general)
     - build

doc/api/index.rst

@@ -227,6 +227,13 @@ and store them in GMT's user data directory.
     datasets.load_earth_vertical_gravity_gradient
     datasets.load_sample_data
 
+In addition, there is also a special function to load XYZ tile maps via
+:doc:`contextily <contextily:index>` to be used as base maps.
+
+.. autosummary::
+    :toctree: generated
+
+    datasets.load_tile_map
 
 .. currentmodule:: pygmt
 

doc/conf.py

@@ -54,11 +54,14 @@
 
 # intersphinx configuration
 intersphinx_mapping = {
-    "python": ("https://docs.python.org/3/", None),
+    "contextily": ("https://contextily.readthedocs.io/en/stable/", None),
     "geopandas": ("https://geopandas.org/en/stable/", None),
     "numpy": ("https://numpy.org/doc/stable/", None),
+    "python": ("https://docs.python.org/3/", None),
     "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None),
+    "rasterio": ("https://rasterio.readthedocs.io/en/stable/", None),
     "xarray": ("https://docs.xarray.dev/en/stable/", None),
+    "xyzservices": ("https://xyzservices.readthedocs.io/en/stable", None),
 }
 
 # options for sphinx-copybutton

doc/install.rst

@@ -106,6 +106,7 @@ PyGMT requires the following libraries to be installed:
 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.
 
 Installing GMT and other dependencies

environment.yml

@@ -12,6 +12,7 @@ dependencies:
     - netCDF4
     - packaging
     # Optional dependencies
+    - contextily
     - geopandas
     - ipython
     # Development dependencies (general)

pygmt/__init__.py

@@ -158,7 +158,15 @@ def _get_ghostscript_version():
         "machine": platform.platform(),
     }
 
-    deps = ["numpy", "pandas", "xarray", "netCDF4", "packaging", "geopandas"]
+    deps = [
+        "numpy",
+        "pandas",
+        "xarray",
+        "netCDF4",
+        "packaging",
+        "contextily",
+        "geopandas",
+    ]
 
     print("PyGMT information:")
     print(f"  version: {__version__}")

pygmt/datasets/__init__.py

@@ -14,3 +14,4 @@
     load_earth_vertical_gravity_gradient,
 )
 from pygmt.datasets.samples import list_sample_data, load_sample_data
+from pygmt.datasets.tile_map import load_tile_map

pygmt/datasets/tile_map.py

@@ -0,0 +1,149 @@
+"""
+Function to load raster tile maps from XYZ tile providers, and load as
+:class:`xarray.DataArray`.
+"""
+
+try:
+    import contextily
+except ImportError:
+    contextily = None
+
+import numpy as np
+import xarray as xr
+
+__doctest_requires__ = {("load_tile_map"): ["contextily"]}
+
+
+def load_tile_map(region, zoom="auto", source=None, lonlat=True, wait=0, max_retries=2):
+    """
+    Load a georeferenced raster tile map from XYZ tile providers.
+
+    The tiles that compose the map are merged and georeferenced into an
+    :class:`xarray.DataArray` image with 3 bands (RGB). Note that the returned
+    image is in a Spherical Mercator (EPSG:3857) coordinate reference system.
+
+    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``) means 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``) means
+        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].
+
+    Returns
+    -------
+    raster : xarray.DataArray
+        Georeferenced 3-D data array of RGB values.
+
+    Raises
+    ------
+    ModuleNotFoundError
+        If ``contextily`` is not installed. Follow
+        :doc:`install instructions for contextily <contextily:index>`, (e.g.
+        via ``pip install contextily``) before using this function.
+
+    Examples
+    --------
+    >>> import contextily
+    >>> from pygmt.datasets import load_tile_map
+    >>> raster = load_tile_map(
+    ...     region=[103.60, 104.06, 1.22, 1.49],  # West, East, South, North
+    ...     source=contextily.providers.Stamen.TerrainBackground,
+    ...     lonlat=True,  # bounding box coordinates are longitude/latitude
+    ... )
+    >>> raster.sizes
+    Frozen({'band': 3, 'y': 1024, 'x': 1536})
+    >>> raster.coords
+    Coordinates:
+      * band     (band) int64 0 1 2
+      * y        (y) float64 1.663e+05 1.663e+05 1.663e+05 ... 1.272e+05 ...
+      * x        (x) float64 1.153e+07 1.153e+07 1.153e+07 ... 1.158e+07 ...
+    """
+    # pylint: disable=too-many-locals
+    if contextily is None:
+        raise ModuleNotFoundError(
+            "Package `contextily` is required to be installed to use this function. "
+            "Please use `pip install contextily` or "
+            "`conda install -c conda-forge contextily` "
+            "to install the package."
+        )
+
+    west, east, south, north = region
+    image, extent = contextily.bounds2img(
+        w=west,
+        s=south,
+        e=east,
+        n=north,
+        zoom=zoom,
+        source=source,
+        ll=lonlat,
+        wait=wait,
+        max_retries=max_retries,
+    )
+
+    # Turn RGBA img from channel-last to channel-first and get 3-band RGB only
+    _image = image.transpose(2, 0, 1)  # Change image from (H, W, C) to (C, H, W)
+    rgb_image = _image[0:3, :, :]  # Get just RGB by dropping RGBA's alpha channel
+
+    # Georeference RGB image into an xarray.DataArray
+    left, right, bottom, top = extent
+    dataarray = xr.DataArray(
+        data=rgb_image,
+        coords={
+            "band": [0, 1, 2],  # Red, Green, Blue
+            "y": np.linspace(start=top, stop=bottom, num=rgb_image.shape[1]),
+            "x": np.linspace(start=left, stop=right, num=rgb_image.shape[2]),
+        },
+        dims=("band", "y", "x"),
+    )
+
+    # If rioxarray is installed, set the coordinate reference system
+    if hasattr(dataarray, "rio"):
+        dataarray = dataarray.rio.write_crs(input_crs="EPSG:3857")
+
+    return dataarray

pyproject.toml

@@ -42,6 +42,7 @@ dynamic = ["version"]
 
 [project.optional-dependencies]
 all = [
+    "contextily",
     "geopandas",
     "ipython"
 ]