Add function to load raster tile maps using contextily

.github/workflows/ci_tests.yaml

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

.github/workflows/ci_tests_dev.yaml

@@ -100,7 +100,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 \
+                        contextily build 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

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

doc/api/index.rst

@@ -227,6 +227,7 @@ and store them in the GMT cache folder.
     datasets.load_fractures_compilation
     datasets.load_hotspots
     datasets.load_japan_quakes
+    datasets.load_map_tiles
     datasets.load_mars_shape
     datasets.load_ocean_ridge_points
     datasets.load_sample_bathymetry

doc/conf.py

@@ -54,11 +54,13 @@
 
 # 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),
     "xarray": ("https://xarray.pydata.org/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
     # Development dependencies (general)
     - build

pygmt/datasets/__init__.py

@@ -4,6 +4,7 @@
 
 from pygmt.datasets.earth_age import load_earth_age
 from pygmt.datasets.earth_relief import load_earth_relief
+from pygmt.datasets.map_tiles import load_map_tiles
 from pygmt.datasets.samples import (
     list_sample_data,
     load_fractures_compilation,

pygmt/datasets/map_tiles.py

@@ -0,0 +1,108 @@
+"""
+Function to load raster basemap tiles 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
+
+
+def load_map_tiles(region, source=None, lonlat=False, **kwargs):
+    """
+    Load a georeferenced raster basemap 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 Spherical Mercator
+        (EPSG:3857) if ``lonlat=False``, or longitude/latitude if
+        ``lonlat=True``.
+
+    source : xyzservices.TileProvider or str
+        [Optional. Default: Stamen Terrain web tiles] The tile source: web tile
+        provider or path to local file. The web tile provider can be in the
+        form of a :class:`xyzservices.TileProvider` object or a URL. The
+        placeholders for the XYZ in the URL need to be {x}, {y}, {z},
+        respectively. For local file paths, the file is read with rasterio and
+        all bands are loaded into the basemap. IMPORTANT: tiles are assumed to
+        be in the Spherical Mercator projection (EPSG:3857).
+
+    lonlat : bool
+        [Optional. Default: False]. If True, coordinates in ``region`` are
+        assumed to be lon/lat as opposed to Spherical Mercator.
+
+    kwargs : dict
+        Extra keyword arguments to pass to :func:`contextily.bounds2img`.
+
+    Returns
+    -------
+    raster : xarray.DataArray
+        Georefenced 3D data array of RGB value.
+
+    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_map_tiles
+    >>> raster = load_map_tiles(
+    ...     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 ...
+    """
+    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, (left, right, bottom, top) = contextily.bounds2img(
+        w=west, s=south, e=east, n=north, source=source, ll=lonlat, **kwargs
+    )
+
+    # Turn RGBA image from channel-last (H, W, C) to channel-first (C, H, W)
+    # and get just RGB (3 band) by dropping RGBA's alpha channel
+    rgb_image = image.transpose(2, 0, 1)[0:3, :, :]
+
+    # Georeference RGB image into an xarray.DataArray
+    dataarray = xr.DataArray(
+        data=rgb_image,
+        coords=dict(
+            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

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