Add Figure.tilemap to plot XYZ tile maps

[weiji14]

Description of proposed changes

New tilemap function for plotting XYZ tile maps! This is essentially a convenience wrapper around pygmt.datasets.load_tile_map and pygmt.Figure.grdimage.

Preview at https://pygmt-dev--2394.org.readthedocs.build/en/2394/api/generated/pygmt.Figure.tilemap.html

This is the second part of #2115 that takes the output xarray.DataArray from pygmt.datasets.load_tile_map (added in #2125) and plots it using fig.grdimage.

Usage:

importcontextilyimportpygmtfig=pygmt.Figure()
fig.tilemap(
region=[-157.6, -157.1, 1.68, 2.08], # West, East, South, Northsource=contextily.providers.Stamen.Watercolor,
frame=True,
)
fig.show()

produces this watercolor map of Kiribati:

Notes:

• rioxarray will be added as an optional dependency, as it is used for the xarray.DataArray to GeoTIFF conversion
• The implementation of tilemap is a little tricky due to cyclic import issues, please suggest if there are better ways fixed with Add Figure.tilemap to plot XYZ tile maps #2394 (comment)
• By default, a given region with lonlat coordinates will produce a map with lonlat (OGC:WGS84) coordinates that are distorted from the original Spherical Mercator (EPSG:3857) tiles, see Add Figure.tilemap to plot XYZ tile maps #2394 (comment)
• The region parameter only accepts a list like [lonmin, lonmax, latmin, latmax], which is different from other GMT/PyGMT functions that accept more types of region values, see Add Figure.tilemap to plot XYZ tile maps #2394 (comment)

Fixes#2115

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If wrapping a new module, open a 'Wrap new GMT module' issue and submit reasonably-sized PRs.
• If adding new functionality, add an example to docstrings or tutorials.
• Use underscores (not hyphens) in names of Python files and directories.

Slash Commands

You can write slash commands (/command) in the first line of a comment to perform
specific operations. Supported slash commands are:

• /format: automatically format and lint the code
• /test-gmt-dev: run full tests on the latest GMT development version

[github-actions]

Summary of changed images

This is an auto-generated report of images that have changed on the DVC remote

Status	Path
added	pygmt/tests/baseline/test_tilemap_no_clip_False.png
added	pygmt/tests/baseline/test_tilemap_no_clip_True.png
added	pygmt/tests/baseline/test_tilemap_ogc_wgs84.png
added	pygmt/tests/baseline/test_tilemap_web_mercator.png

Image diff(s)

Details

Added images

• pygmt/tests/baseline/test_tilemap_no_clip_False.png

• pygmt/tests/baseline/test_tilemap_no_clip_True.png

• pygmt/tests/baseline/test_tilemap_ogc_wgs84.png

• pygmt/tests/baseline/test_tilemap_web_mercator.png

Modified images

Path	Old	New

Report last updated at commit c3fca14

[weiji14]

Currently region needs to be a list to work with load_tile_map, but I'm debating on whether to allow strings like 0/1/2/3 or other styles supported by GMT -R as mentioned in #153 (comment) and https://docs.generic-mapping-tools.org/6.4/cookbook/options.html#data-domain-or-map-region-the-r-option

[seisman]

other styles supported by GMT -R as mentioned in #153 (comment) and https://docs.generic-mapping-tools.org/6.4/cookbook/options.html#data-domain-or-map-region-the-r-option

Other styles are more tricky to support. For example, how can we know the boundingbox for region="US"?

[weiji14]

Other styles are more tricky to support.

Hmm, I was mostly just thinking about 0/1/2/3, and maybe np.array([0, 1, 2, 3]). Usually users would set a region variable near the top of their code, and they might find it strange why their region that works for fig.basemap doesn't work for fig.tilemap.

For example, how can we know the boundingbox for region="US"?

If only there was an API in GMT to return the bounding box for a given region 😉 Global regions like region="g" and region="d" might be easy to support, but ISO code regions (https://www.pygmt.org/v0.8.0/tutorials/basics/regions.html#iso-code) like region="JP" and `region=JP+r3" will be a bit trickier.

[seisman]

If only there was an API in GMT to return the bounding box for a given region Global regions like region="g" and region="d" might be easy to support, but ISO code regions (https://www.pygmt.org/v0.8.0/tutorials/basics/regions.html#iso-code) like region="JP" and `region=JP+r3" will be a bit trickier.

Open a upstream feature request in the GMT repository?

[weiji14]

Hmm, or do you think we can use GMT_Extract_Region that's wrapped here:

pygmt/pygmt/clib/session.py

Lines 1500 to 1571 in 0b5a69b

	defextract_region(self):
	"""
	Extract the WESN bounding box of the currently active figure.
	Retrieves the information from the PostScript file, so it works for
	country codes as well.
	Returns
	-------
	* wesn : 1-D array
	A numpy 1-D array with the west, east, south, and north dimensions
	of the current figure.
	Examples
	--------
	>>> import pygmt
	>>> fig = pygmt.Figure()
	>>> fig.coast(
	... region=[0, 10, -20, -10],
	... projection="M6i",
	... frame=True,
	... land="black",
	... )
	>>> with Session() as lib:
	... wesn = lib.extract_region()
	...
	>>> print(", ".join([f"{x:.2f}" for x in wesn]))
	0.00, 10.00, -20.00, -10.00
	Using ISO country codes for the regions (for example ``'US.HI'`` for
	Hawaii):
	>>> fig = pygmt.Figure()
	>>> fig.coast(
	... region="US.HI", projection="M6i", frame=True, land="black"
	... )
	>>> with Session() as lib:
	... wesn = lib.extract_region()
	...
	>>> print(", ".join([f"{x:.2f}" for x in wesn]))
	-164.71, -154.81, 18.91, 23.58
	The country codes can have an extra argument that rounds the region a
	multiple of the argument (for example, ``'US.HI+r5'`` will round the
	region to multiples of 5):
	>>> fig = pygmt.Figure()
	>>> fig.coast(
	... region="US.HI+r5", projection="M6i", frame=True, land="black"
	... )
	>>> with Session() as lib:
	... wesn = lib.extract_region()
	...
	>>> print(", ".join([f"{x:.2f}" for x in wesn]))
	-165.00, -150.00, 15.00, 25.00
	"""
	c_extract_region=self.get_libgmt_func(
	"GMT_Extract_Region",
	argtypes=[ctp.c_void_p, ctp.c_char_p, ctp.POINTER(ctp.c_double)],
	restype=ctp.c_int,
	)
	wesn=np.empty(4, dtype=np.float64)
	wesn_pointer=wesn.ctypes.data_as(ctp.POINTER(ctp.c_double))
	# The second argument to GMT_Extract_Region is a file pointer to a
	# PostScript file. It's only valid in classic mode. Use None to get a
	# NULL pointer instead.
	status=c_extract_region(self.session_pointer, None, wesn_pointer)
	ifstatus!=0:
	raiseGMTCLibError("Failed to extract region from current figure.")
	returnwesn

Downside is that it requires a figure instance to be set up already first, before the bounding box coordinates can be extracted.

[seisman]

I don't think it will work. Actually, a figure instance is not enough. We have to plot something to generate the hidden PostScript file before calling extract_region.

In [1]: import pygmt
In [2]: fig = pygmt.Figure()
In [3]: fig.region
pygmt-session [ERROR]: No hidden PS file found
---------------------------------------------------------------------------
GMTCLibError Traceback (most recent call last)
Cell In[3], line 1
----> 1 fig.region
File ~/OSS/gmt/pygmt/pygmt/figure.py:126, in Figure.region(self)
124 self._activate_figure()
125 with Session() as lib:
--> 126 wesn = lib.extract_region()
127 return wesn
File ~/OSS/gmt/pygmt/pygmt/clib/session.py:1570, in Session.extract_region(self)
1568 status = c_extract_region(self.session_pointer, None, wesn_pointer)
1569 if status != 0:
-> 1570 raise GMTCLibError("Failed to extract region from current figure.")
1571 return wesn

[weiji14]

I don't think it will work. Actually, a figure instance is not enough. We have to plot something to generate the hidden PostScript file before calling extract_region.

Hmm yeah, we would need to maybe do a fig.basemap and then fig.grdimage, but that's a bit of a hacky solution.

Let's just keep this PR simple and support only a list of [xmin, xmax, ymin, ymax] coordinates for now. There can be a follow up PR to support other region arguments (most likely after PyGMT v0.9.0) if there's interest in this 'enhancement'.

[weiji14]

Note that this is not the full implementation of PyGMT's fig.grdimage, which is actually a bit more complicated to support shading (-I) with an xarray.DataArray grid

pygmt/pygmt/src/grdimage.py

Lines 179 to 192 in 8f31706

	withSession() aslib:
	file_context=lib.virtualfile_from_data(check_kind="raster", data=grid)
	withcontextlib.ExitStack() asstack:
	# shading using an xr.DataArray
	ifkwargs.get("I") isnotNoneanddata_kind(kwargs["I"]) =="grid":
	shading_context=lib.virtualfile_from_data(
	check_kind="raster", data=kwargs["I"]
	)
	kwargs["I"] =stack.enter_context(shading_context)
	fname=stack.enter_context(file_context)
	lib.call_module(
	module="grdimage", args=build_arg_string(kwargs, infile=fname)
	)

Question is: should we support shading of the XYZ tiles with xarray.DataArray grids? Or is there a way to use the grdimage implementation above directly without code duplication?

[michaelgrund]

Currently I would say no to shading.

[seisman]

other styles supported by GMT -R as mentioned in #153 (comment) and https://docs.generic-mapping-tools.org/6.4/cookbook/options.html#data-domain-or-map-region-the-r-option

Other styles are more tricky to support. For example, how can we know the boundingbox for region="US"?

[weiji14]

Was playing with the zoom parameter just now, and it seems that setting lower zoom levels (e.g. 0) will actually cause the area covered by contextily.bounds2img to be bigger than the original bounding box region. This might come as a surprise for users, so maybe we should set no_clip to False by default? Would require a bit of extra work to get correct. Edit: Sorry, got a bit confused by the double negatives, no_clip (-N) is already False by default, just needed to pass the region arguments to grdimage properly.

Alternatively, we could also do the clipping in pygmt.datasets.load_tile_map.

[weiji14]

Added commit 17d0734 to ensure that:

• The plot is clipped to the bounding box region when no_clip=False (default)
• The plot extends to areas outside of the bounding box region when no_clip=True

[seisman]

Usage:

importcontextilyimportpygmtfig=pygmt.Figure()
fig.tilemap(
region=[-157.6, -157.1, 1.68, 2.08], # West, East, South, Northsource=contextily.providers.Stamen.Watercolor,
frame=True,
)
fig.show()

produces this watercolor map of Kiribati:

Do we need to convert the map tiles from EPSG:3857 to longitude/latitude?

[weiji14]

Do we need to convert the map tiles from EPSG:3857 to longitude/latitude?

Yes, if lonlat=True (default), then I think the output map should be in lon/lat as well to make it intuitive. I've implemented this in fde94ad. The issue with this is that the pixels will be reprojected and distorted by default, which might not be a problem for tile maps like Stamen Terrain, but can be problematic when there are words on the tile map (e.g. for OpenStreetMap).

[seisman]

Setting raster.gmt.gtype = 1 is useless, because raster.rio.to_raster doesn't understand it, right?

[weiji14]

Yes, rioxarray's .rio.to_raster won't understand .gmt.gtype, but I don't want to forget to add this after #2398 is merged 🙂

[seisman]

When no_clip is not given, it defaults to None and not kwargs.get("N") is True, which means the plot is clipped. I find the code is a little difficult to understand. Do you think the following changes make it better? Or if kwargs.get("N") is in [None, False]?

Suggested change

	# Only set region if no_clip is False, so that plot is clipped to exact
	# bounding box region
	ifnotkwargs.get("N"):
	kwargs["R"] ="/".join(str(coordinate) forcoordinateinregion)
	# Only set region if no_clip is False or None, so that plot is clipped to exact
	# bounding box region
	ifkwargs.get("N")isnotTrue:
	kwargs["R"] ="/".join(str(coordinate) forcoordinateinregion)

[weiji14]

Yeah, double negatives are confusing. I went with kwargs.get("N") in [None, False] at 125574e.

[seisman]

This PR looks great to me.

[seisman]

@michaelgrund and @yvonnefroehlich It would be great if you can give this PR a final review before we merge it.

[michaelgrund]

Looks good! Thanks again for working on that @weiji14 😉

[yvonnefroehlich]

I am afriad I do not have the knowledge to review this PR reliable, but what I understand make sense to me.

[seisman]

@weiji14 OK to merge?

[weiji14]

@weiji14 OK to merge?

Yes I'll merge it today. Let me incorporate #2394 (comment) first.

[weiji14]

Hold on, let me fix this dtype error on Windows / Python 3.11 + NumPy 1.24 at https://github.com/GenericMappingTools/pygmt/actions/runs/4527632070/jobs/7973665357?pr=2394#step:11:737. Looks similar to #2393:

_______________ [doctest] pygmt.datasets.tile_map.load_tile_map _______________
095 >>> from pygmt.datasets import load_tile_map
096 >>> raster = load_tile_map(
097 ... region=[-180.0, 180.0, -90.0, 0.0], # West, East, South, North
098 ... zoom=1, # less detailed zoom level
099 ... source=contextily.providers.Stamen.TerrainBackground,
100 ... lonlat=True, # bounding box coordinates are longitude/latitude
101 ... )
102 >>> raster.sizes
103 Frozen({'band': 3, 'y': 256, 'x': 512})
104 >>> raster.coords
Differences (unified diff with -expected +actual):
@@-1,5+1,5@@
Coordinates:
* band (band) uint8 012-* y (y) float64 -7.081e-10-7.858e+04...-1.996e+07...+* 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+072.004e+07- spatial_ref int64 0+ spatial_ref int32 0
D:\a\pygmt\pygmt\pygmt\datasets\tile_map.py:104: DocTestFailure
============================== warnings summary ===============================

[weiji14]

Changed from using .rio.write_crs to .rio.set_crs to prevent the CF coordinate spatial_ref from being written to the output xarray.DataArray. This fixes the int32 on Windows bu int64 on Linux/macOS issue mentioned in #2394 (comment).

If this looks ok (and the tests pass), I'll proceed to merge in this PR.

[seisman]

Looks good.