diff --git a/.github/workflows/ci_docs.yml b/.github/workflows/ci_docs.yml index 5ec85256512..c7a87a005bc 100644 --- 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 index 974f46a3286..6cd6b0a1c5b 100644 --- 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 index d685d261a50..bbeaba52c44 100644 --- 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 index 2b2d928b06c..fd47dae75d4 100644 --- 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 index 7d725cde8e2..a11194fec36 100644 --- 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 index cf017a475bd..6303559be52 100644 --- 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 ` 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 index fd59db51c9a..9fa464c5c67 100644 --- 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 index 09c7aa9d968..87efd1039d5 100644 --- 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 `__: For embedding the figures in Jupyter notebooks (recommended). +* `Contextily `__: For retrieving tile maps from the internet. * `GeoPandas `__: For using and plotting GeoDataFrame objects. Installing GMT and other dependencies diff --git a/environment.yml b/environment.yml index fcbfea58e47..6f7ca48edb4 100644 --- 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 index deec3f2b114..1e1ee13447b 100644 --- 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 index d7a099972fb..7aba16154be 100644 --- 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 new file mode 100644 index 00000000000..2afe28acbf9 --- /dev/null +++ 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 ` 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 ` 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 `, (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 index ac52fba797a..f0880e22304 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,6 +42,7 @@ dynamic = ["version"] [project.optional-dependencies] all = [ + "contextily", "geopandas", "ipython" ]