Integration with contextily to plot xyz basemaps

[weiji14]

Description of the desired feature

To allow for a greater variety of basemaps from e.g. OpenStreetMap, Stamen, MapBox, etc, it would be nice to integrate with contextily (part of the geopandas organization). This is motivated by https://forum.generic-mapping-tools.org/t/pygmt-vs-plotly/3250/10.

On the implementation, here are 2 ideas on how to integrate contextily basemaps with pygmt.

1. Have an additional option under fig.basemap which allows selection of the basemap (e.g. OpenStreetMap, Stamen, Mapbox, etc)
2. Create a new load_basemap function under pygmt/datasets that would load the image (np.ndarray) and georeference it (to an xarray.DataArray), and then the user plots it using fig.grdimage.
3. Maybe there's another more Pythonic way? Comment down below!

Either way, these would require:

1. Wrapping around the contextily.bounds2img function at https://contextily.readthedocs.io/en/latest/reference.html#contextily.bounds2img. This returns a numpy.ndarray and a geographic extent.
2. The numpy.ndarray needs to be georeferenced into an xarray.DataArray. I have some recent code at https://zen3geo.readthedocs.io/en/v0.4.0/object-detection-boxes.html#georeference-image-using-rioxarray to do that
3. The xarray.DataArray will need to be plotted on the PyGMT map

Are you willing to help implement and maintain this feature? Yes, but discuss about the implementation first

[seisman]

On the implementation, here are 2 ideas on how to integrate contextily basemaps with pygmt.

1. Have an additional option under fig.basemap which allows selection of the basemap (e.g. OpenStreetMap, Stamen, Mapbox, etc)
2. Create a new load_basemap function under pygmt/datasets that would load the image (np.ndarray) and georeference it (to an xarray.DataArray), and then the user plots it using fig.grdimage.
3. Maybe there's another more Pythonic way? Comment down below!

Can we do both 1 and 2? I.e., having the load_basemap (maybe load_map_tiles?) function in pygmt/datasets and having a new parameter in Figure.basemap to call this function.

[weiji14]

On the implementation, here are 2 ideas on how to integrate contextily basemaps with pygmt.

1. Have an additional option under fig.basemap which allows selection of the basemap (e.g. OpenStreetMap, Stamen, Mapbox, etc)
2. Create a new load_basemap function under pygmt/datasets that would load the image (np.ndarray) and georeference it (to an xarray.DataArray), and then the user plots it using fig.grdimage.
3. Maybe there's another more Pythonic way? Comment down below!

Can we do both 1 and 2? I.e., having the load_basemap (maybe load_map_tiles?) function in pygmt/datasets and having a new parameter in Figure.basemap to call this function.

That sounds like a neat idea. Start with (2) load_map_tiles first, and then do (1) new option in fig.basemap() (which will be a lot more complicated.

Actually, we'll also need to let fig.grdimage work with RGB (3 bands) instead of just 1 band. It works if we save the RGB image to a GeoTIFF first and tell fig.grdimage to plot the .tif file, but would be nicer to plot a 3 band xr.DataArray directly, xref #1555. Or we can use a tempfile 😉

[weiji14]

Started a Pull Request at #2125 for ~load_map_tiles~ load_tile_map (renamed based on https://github.com/GenericMappingTools/pygmt/pull/2125#discussion_r1089730111).

To jump ahead a little bit, this is the code I'm currently using to plot the 3-band xarray.DataArray. The raster can be obtained from running the snippet in https://github.com/GenericMappingTools/pygmt/pull/2125#issue-1378213857.

import pygmt
import rioxarray

raster.rio.to_raster(raster_path="basemap.tif")

fig = pygmt.Figure()
fig.grdimage(grid="basemap.tif", frame=True)
fig.savefig("Singapore_stamen_map.png")
fig.show()

produces this Stamen Terrain map over Singapore

[weiji14]

1. Have an additional option under fig.basemap which allows selection of the basemap (e.g. OpenStreetMap, Stamen, Mapbox, etc)

After a bit of thinking, and as mentioned in in https://github.com/GenericMappingTools/pygmt/pull/2125#discussion_r1059773318, I'm considering creating a new fig.tilemap method instead of complicating fig.basemap further. This fig.tilemap method would essentially look like:

def tilemap(self, *, **kwargs):
    with GMTTempFile(suffix=".tif") as tmpfile:
        raster = pygmt.datasets.load_tile_map(...)
        raster.rio.to_raster(raster_path=tmpfile.name)
        fig.grdimage(grid=tmpfile.name, **kwargs)

Any thoughts on this implementation? Note that projections will also need to be handled, something mentioned in https://github.com/GenericMappingTools/pygmt/pull/2125#issuecomment-1447292116. Also, the .rio.to_raster requires adding rioxarray as another optional dependency.

The tempfile is a bit of a hack to be honest, but is necessary until we can use grdimage to plot 3-band xarray.DataArrays directly, something to be done in #1555.

[seisman]

After a bit of thinking, and as mentioned in in #2125 (comment), I'm considering creating a new fig.tilemap method instead of complicating fig.basemap further.

Yes, I aslo prefer fig.tilemap!

The tempfile is a bit of a hack to be honest, but is necessary until we can use grdimage to plot 3-band xarray.DataArrays directly, something to be done in #1555.

I agree.