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