load_tile_map: Add the new parameter 'crs' to set the CRS of the returned dataarray

[seisman]

Description of proposed changes

This PR adds the crs parameter to the load_tile_map function so that users can change the CRS of the returned raster image without needing to know the details.

Address #2125 (comment).
Closes #3484.

Notes for maintainers:

As far as I know, there are multiple different ways to reproject raster images.

1. contextily.warp_tiles: Internally calls rasterio.vrt.WarpedVRT
2. rioxarray's rio.reproject: Internally calls rasterio.warp.reproject
3. rasterio.vrt.WarpedVRT
4. rasterio.warp.reproject
5. gdal.Warp
6. ...

Options 3-5 are complicated and I don't think we want to call them directly. Ideally, we should use option 1, so reprojection doesn't rely on an extra dependency rioxarray. However, as shown below, the result from contextily.warp_tiles doesn't work well with rio.to_raster.

Here are codes to compare the results with options 1 and 2:

import xarray as xr
import rioxarray
import contextily as cx
import numpy as np

def image2raster(image, extent, crs="EPSG:3857"):
    """
    Convert an image to a xarray.DataArray raster.

    Copied from pygmt.datasets.tile_map.load_tile_map, except that crs is an input parameter.
    """
    # 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": np.array(object=[1, 2, 3], dtype=np.uint8),  # 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"),
    )
    dataarray = dataarray.rio.write_crs(input_crs=crs)
    return dataarray

# Image and extent of the global tile
image, extent = cx.bounds2img(w=-180, e=180, s=-90, n=90, ll=True, zoom=0)

# Option 1: Use contextily.warp_tiles
image_warp, extent_warp = cx.warp_tiles(image, extent, "OGC:CRS84")
dataarray_warp = image2raster(image_warp, extent_warp, crs="OGC:CRS84")
print("Using contextily.warp_tiles:", dataarray_warp.rio.bounds())
dataarray_warp.rio.to_raster("result_warp.tiff")

# Option 2: Use rio.reproject
dataarray = image2raster(image, extent, crs="EPSG:3857")
dataarray_reproject = dataarray.rio.reproject("OGC:CRS84")
print("Using rio.reproject:", dataarray_reproject.rio.bounds())
dataarray_reproject.rio.to_raster("result_reproject.tiff")

The outputs are:

Using contextily.warp_tiles: (-180.55157789333614, -85.9082903958547, 180.18036434850873, 85.66487761291295)
Using rio.reproject: (-180.0, -84.8307625666717, 180.00000000000006, 85.11165075221737)

As you can see, the longitude range is (-180.55157789333614, 180.18036434850873) with contextily.warp_tiles, which will lead to GMT errors like:

grdimage [ERROR]: Map region exceeds 360 degrees
grdimage [ERROR]: General map projection error

Two more notes:

• There is an open request to use rasterio.warp.reproject instead of rasterio.vrt.WarpedVRT in warp_tiles (Explore rasterio.warp.reproject to replace MemoryFile for warping geopandas/contextily#88). If it's done, then the result from warp_tiles will be likely the same as the one from rio.reproject
• Currently, we need to save raster into temporary files by calling rio.to_raster (

  pygmt/pygmt/helpers/tempfile.py

  Line 204 in 988746b

  image.rio.to_raster(raster_path=tmpfile.name)

  ). We don't have to do it after WIP: Implement virtualfile_from_image #3468.

Preview: https://pygmt-dev--3554.org.readthedocs.build/en/3554/api/generated/pygmt.datasets.load_tile_map.html

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 command is:

• /format: automatically format and lint the code

[seisman]

The Figure.tilemap method no longer requires rioxarray explicitly. Instead, load_tile_map will raise the error.

[seisman]

If lonlat=True, set the CRS to let load_tile_map do the raster reprojection for us.

[michaelgrund]

Looks good to me!

[weiji14]

Maybe put the rasterio imports in the rioxarray block below instead of here.

[seisman]

rasterio is a known dependency of contextily (xref: https://github.com/geopandas/contextily/blob/main/pyproject.toml), so it makes more sense to assume that they're installed together.

[weiji14]

Rasterio is also a dependency of rioxarray 😆. I would argue that rasterio is more tightly coupled to rioxarray, and since rasterio/rioxarray are both used for their projection system capabilities, it would make more sense to put those together. But ok too if you want to keep it here.

[weiji14]

At https://contextily.readthedocs.io/en/latest/reference.html#contextily.bounds2img, there is this sentence:

IMPORTANT: tiles are assumed to be in the Spherical Mercator projection (EPSG:3857), unless the crs keyword is specified.

While tilemaps are usually served in EPSG:3857, I know that there are some tilemaps that can be served in a different projection system by default. I'm hesitant to use this default of EPSG:3857 here, unless we can be confident that contextily only supports EPSG:3857 tiles.

[seisman]

Looking at the contexily source codes, I think contexitly always assumes that CRS is EPSG:3857, at least when writing the image to raster files

https://github.com/geopandas/contextily/blob/eafeb3c993ab43356b65e7a29b9dd6df3ac76f00/contextily/tile.py#L358

[weiji14]

That is for the warp_tiles function. I don't see anything about EPSG:3857 in the bounds2img function (though it's calling several functions and I haven't checked too closely).

Non-Web Mercator (EPSG:3857) XYZ tiles are not common, but they do exist, e.g. Openlayers allows configuring the projection system of the returned tiles:

• https://www.linz.govt.nz/guidance/data-service/linz-data-service-guide/map-tile-services/using-lds-xyz-services-openlayers#c-0-s-4
• https://openlayers.org/en/main/apidoc/module-ol_source_XYZ-XYZ.html
• https://stackoverflow.com/questions/65024840/xyz-in-non-web-mercator-projections

[seisman]

Looking at the contexily source codes, I think contexitly always assumes that CRS is EPSG:3857, at least when writing the image to raster files

I meant this line:

https://github.com/geopandas/contextily/blob/eafeb3c993ab43356b65e7a29b9dd6df3ac76f00/contextily/tile.py#L175

In bounds2raster, it first calls bounds2img to get the image and extent (can be any CRS), but when writing to a raster image, bounds2raster assumes the image is EPSG:3857.

Looking at the xyzservices (https://github.com/geopandas/xyzservices/blob/c847f312d2e67a0a45ceae119fdcdf7cf60858b6/provider_sources/xyzservices-providers.json#L50), I can also see some providers that has the crs attribute, but it's unclear how the crs attribute is handled in contexitly.

[weiji14]

I tried using xyzservices.providers.MapTiler.Basic4326 which comes in EPSG:4326:

import contextily
import xyzservices
import matplotlib.pyplot as plt

source = xyzservices.providers.MapTiler.Basic4326(
    key="BwTZdryxoypHc4BnCZ4"  # add an 'o' to the end
)
img, bounds = contextily.bounds2img(w=0, s=-80, e=150, n=80, ll=True, source=source)
# img.shape # (4096, 2048, 4)
plt.imshow(img[:, :, :3])

produces

The extent doesn't seem to be correct (it should be plotting from Greenwich to Asia), possibly a bug in contextily hardcoding to EPSG:3857. But it does show that non-Web Mercator XYZ tiles are indeed possible.

[weiji14]

Found geopandas/contextily#119, and had a closer look at the bounds2img code, specifically these lines - https://github.com/geopandas/contextily/blob/00ca0318033e69e29a164ce571d9c932a057ecc6/contextily/tile.py#L270-L272. It seems like contextily is using mercantile which only supports Spherical Mercator (EPSG:3857), see mapbox/mercantile#109 (comment). Unless contextily switches to something like https://github.com/developmentseed/morecantile in the future, then it would only 'work' with EPSG:3857.

[seisman]

The extent doesn't seem to be correct (it should be plotting from Greenwich to Asia), possibly a bug in contextily hardcoding to EPSG:3857.

https://github.com/geopandas/contextily/blob/eafeb3c993ab43356b65e7a29b9dd6df3ac76f00/contextily/tile.py#L256

When ll=True, bounds2img calls the private _sm2ll function to convert bounds in Spherical Mercator coordinates to lon/lat. I believe "EPSG:3857" is hardcoded/assumed in many places in the contextily package.

[weiji14]

How about checking source.crs (an attribute from an instance of the xyzservices.TileProvider class), and if it is present, use that as the default crs. If not, fallback to EPSG:3857.

[seisman]

Sounds good. Done in 6733e19.

_default_crs = getattr(source, "crs", "EPSG:3857")

Please note that this line of code works for TileProvide, str and None. str and None don't have crs, so EPSG:3857 is returned.

Edit: renamed _default_crs to _source_crs in c1d55b0.

Please note that, currently, the returned xarray.DataArray is still "EPSG:3857" even if the source CRS is not (i.e. we don't the reprojection). This may not be ideal.

[weiji14]

That is for the warp_tiles function. I don't see anything about EPSG:3857 in the bounds2img function (though it's calling several functions and I haven't checked too closely).

Non-Web Mercator (EPSG:3857) XYZ tiles are not common, but they do exist, e.g. Openlayers allows configuring the projection system of the returned tiles:

• https://www.linz.govt.nz/guidance/data-service/linz-data-service-guide/map-tile-services/using-lds-xyz-services-openlayers#c-0-s-4
• https://openlayers.org/en/main/apidoc/module-ol_source_XYZ-XYZ.html
• https://stackoverflow.com/questions/65024840/xyz-in-non-web-mercator-projections