diff --git a/.github/workflows/ci_docs.yml b/.github/workflows/ci_docs.yml index ad6587000a0..9a3306a620a 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 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 diff --git a/.github/workflows/ci_tests.yaml b/.github/workflows/ci_tests.yaml index 8a9b5d94c8d..7e0ca7995b7 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: 'contextily geopandas ipython' + optional-packages: 'contextily geopandas ipython rioxarray' timeout-minutes: 30 defaults: run: diff --git a/.github/workflows/ci_tests_dev.yaml b/.github/workflows/ci_tests_dev.yaml index fa73d1c0c4f..c5b19178d41 100644 --- a/.github/workflows/ci_tests_dev.yaml +++ b/.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 diff --git a/.github/workflows/ci_tests_legacy.yaml b/.github/workflows/ci_tests_legacy.yaml index c956e3a27fb..5a1506327fd 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 \ - contextily geopandas ipython \ + contextily geopandas ipython rioxarray \ 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 31decb7123f..f623b310bc1 100644 --- a/ci/requirements/docs.yml +++ b/ci/requirements/docs.yml @@ -14,6 +14,7 @@ dependencies: # Optional dependencies - contextily - geopandas + - rioxarray # Development dependencies (general) - build - ipython diff --git a/doc/api/index.rst b/doc/api/index.rst index 8bd68a01d42..29c1ea3bdf1 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -62,6 +62,7 @@ Plotting raster data Figure.grdimage Figure.grdview Figure.image + Figure.tilemap Configuring layout ~~~~~~~~~~~~~~~~~~ diff --git a/doc/conf.py b/doc/conf.py index 9fa464c5c67..6c356c0fe2e 100644 --- a/doc/conf.py +++ b/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), } diff --git a/doc/install.rst b/doc/install.rst index 9d45ca5ff5b..834fb168ae3 100644 --- a/doc/install.rst +++ b/doc/install.rst @@ -108,6 +108,7 @@ 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. +* `RioXarray `__: For saving multi-band rasters to GeoTIFFs. Installing GMT and other dependencies ------------------------------------- diff --git a/environment.yml b/environment.yml index 3ed7b19fd2d..3ddd497e019 100644 --- a/environment.yml +++ b/environment.yml @@ -15,6 +15,7 @@ dependencies: - contextily - geopandas - ipython + - rioxarray # Development dependencies (general) - build - dvc diff --git a/pygmt/datasets/tile_map.py b/pygmt/datasets/tile_map.py index 260232ba197..fd5268a408a 100644 --- a/pygmt/datasets/tile_map.py +++ b/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 diff --git a/pygmt/figure.py b/pygmt/figure.py index 728fbf3a886..b18d5bac974 100644 --- a/pygmt/figure.py +++ b/pygmt/figure.py @@ -523,6 +523,7 @@ def _repr_html_(self): subplot, ternary, text, + tilemap, timestamp, velo, wiggle, diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py index 087e8a063b4..e2fd1752750 100644 --- a/pygmt/src/__init__.py +++ b/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 diff --git a/pygmt/src/tilemap.py b/pygmt/src/tilemap.py new file mode 100644 index 00000000000..be2338fa05e --- /dev/null +++ b/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 ` 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``]. + + 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 `, + (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) + ) diff --git a/pygmt/tests/baseline/test_tilemap_no_clip_False.png.dvc b/pygmt/tests/baseline/test_tilemap_no_clip_False.png.dvc new file mode 100644 index 00000000000..8de7c7e7197 --- /dev/null +++ b/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 diff --git a/pygmt/tests/baseline/test_tilemap_no_clip_True.png.dvc b/pygmt/tests/baseline/test_tilemap_no_clip_True.png.dvc new file mode 100644 index 00000000000..e5c1f4935ea --- /dev/null +++ b/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 diff --git a/pygmt/tests/baseline/test_tilemap_ogc_wgs84.png.dvc b/pygmt/tests/baseline/test_tilemap_ogc_wgs84.png.dvc new file mode 100644 index 00000000000..94a2d317aa3 --- /dev/null +++ b/pygmt/tests/baseline/test_tilemap_ogc_wgs84.png.dvc @@ -0,0 +1,4 @@ +outs: +- md5: 3de0555d86aca49b92425c8d5272a934 + size: 59286 + path: test_tilemap_ogc_wgs84.png diff --git a/pygmt/tests/baseline/test_tilemap_web_mercator.png.dvc b/pygmt/tests/baseline/test_tilemap_web_mercator.png.dvc new file mode 100644 index 00000000000..ad6cc210818 --- /dev/null +++ b/pygmt/tests/baseline/test_tilemap_web_mercator.png.dvc @@ -0,0 +1,4 @@ +outs: +- md5: a76d9a9a1890d6b1345305eaea598bc3 + size: 122195 + path: test_tilemap_web_mercator.png diff --git a/pygmt/tests/test_tilemap.py b/pygmt/tests/test_tilemap.py new file mode 100644 index 00000000000..06d377f4462 --- /dev/null +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index ac3c37ad57d..bdd62c8e28b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ keywords = [ "geophysics", "geospatial", "oceanography", - "seismology" + "seismology", ] classifiers = [ "Development Status :: 4 - Beta", @@ -36,7 +36,7 @@ dependencies = [ "pandas", "xarray", "netCDF4", - "packaging" + "packaging", ] dynamic = ["version"] @@ -44,7 +44,8 @@ dynamic = ["version"] all = [ "contextily", "geopandas", - "ipython" + "ipython", + "rioxarray", ] [project.urls]