Add function to load raster tile maps using contextily (#2125)

weiji14 · seisman · michaelgrund · web-flow · commit f379164d1b3f · 2023-03-03T11:31:23.000+13:00
New dataset function to load XYZ tiles! Uses contextily to retrieve the tiles based on a bounding box. Included an example doctest that shows how the map tiles can be loaded into an xarray.DataArray. Added a new section in the API docs and intersphinx mappings for contextily and xyzservices. * Use correct Spherical Mercator coordinates Can't assume that the input bounding box (which can be in longitude/latitude) is the same as the returned extent (which is always in EPSG:3857). * Change ll parameter to lonlat To fix pylint `C0103: Argument name "ll" doesn't conform to snake_case naming style (invalid-name)`. * Add contextily to CI build matrix and include it as optional dependency Let the Continuous Integration tests run with `contextily`, include it in pyproject.toml and environment.yml, and document it in `doc/install.rst` as an optional dependency. * Set default bounding box coordinates to be lonlat Bounding box coordinates are assumed to be longitude/latitude by default, rather than in Spherical Mercator. * Skip doctest when contextily is not installed Using the `__doctest_requires__` variable, see https://github.com/astropy/pytest-doctestplus/tree/v0.12.1#doctest-dependencies * Add intersphinx link for rasterio Also updated intersphinx link for xarray to new URL at https://docs.xarray.dev/en/stable * Document wait and max_retries parameters used in contextily.bounds2img * Use PyGMT's convention for default values in docstrings Modified from original contextily.bounds2img docstring to fit PyGMT's standards. Xref https://github.com/Generi cMappingTools/pygmt/pull/1182. * Rename load_map_tiles to load_tile_map Also create new dedicated section for load_tile_map in the API docs index. * Add zoom parameter and remove kwargs Wrap all of the parameters in contextily.bounds2img, and so can remove kwargs. Need to disable the pylint recommendation `R0914: Too many local variables`. * Add contextily to docs build CI requirements * Add contextily to pygmt.show_versions() dependency list * Document the three possible source options thoroughly Split out the three possible source options into bullet points. Link to https://contextily.readthedocs.io/en/latest/providers_deepdive.html, give an example OpenStreetMap URL, and link to https://contextily.readthedocs.io/en/stable/working_with_local_files.html. * Add more detail about the zoom level of detail --------- Co-authored-by: Dongdong Tian <seisman.info@gmail.com> Co-authored-by: Michael Grund <23025878+michaelgrund@users.noreply.github.com> Co-authored-by: Yvonne Fröhlich <94163266+yvonnefroehlich@users.noreply.github.com>
diff --git a/.github/workflows/ci_docs.yml b/.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
diff --git a/.github/workflows/ci_tests.yaml b/.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:
diff --git a/.github/workflows/ci_tests_dev.yaml b/.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)
diff --git a/.github/workflows/ci_tests_legacy.yaml b/.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
 
diff --git a/ci/requirements/docs.yml b/ci/requirements/docs.yml
@@ -12,6 +12,7 @@ dependencies:
     - netCDF4
     - packaging
     # Optional dependencies
+    - contextily
     - geopandas
     # Development dependencies (general)
     - build
diff --git a/doc/api/index.rst b/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
 
diff --git a/doc/conf.py b/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
diff --git a/doc/install.rst b/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
diff --git a/environment.yml b/environment.yml
@@ -12,6 +12,7 @@ dependencies:
     - netCDF4
     - packaging
     # Optional dependencies
+    - contextily
     - geopandas
     - ipython
     # Development dependencies (general)
diff --git a/pygmt/__init__.py b/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__}")
diff --git a/pygmt/datasets/__init__.py b/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
diff --git a/pygmt/datasets/tile_map.py b/pygmt/datasets/tile_map.py
@@ -0,0 +1,151 @@
+"""
+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``) 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``].
+
+    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
diff --git a/pyproject.toml b/pyproject.toml
@@ -42,6 +42,7 @@ dynamic = ["version"]
 
 [project.optional-dependencies]
 all = [
+    "contextily",
     "geopandas",
     "ipython"
 ]
