Figure.tilemap: handle longitude given as 0°-360°

[yvonnefroehlich]

pygmt.Figure.coast and the remote datasets (plotted via pygmt.Figure.grdimage) support that the longitude is given either as -180° to 180° East or as 0° to 360° East. However, it looks like the later longitude range is not supported by pygmt.Figure.tilemap. When plotting tiles for a region exceeding 180° East, the area after 180° East remains empty.

import pygmt as gmt

region_lon360 = [120, 190, -48, -9]

fig = gmt.Figure()

fig.tilemap(
    region=region_lon360,
    projection="M10c",
    zoom=2,
    frame=["WSne+ttilemap", "af"],
)

fig.shift_origin(xshift="11c")

fig.coast(
    region=region_lon360,
    projection="M10c",
    shorelines="1/1p,black",
    water="steelblue",
    frame=["wSne+tcoast", "af"],
)

fig.shift_origin(xshift="11c")

fig.grdimage(
    region=region_lon360,
    projection="M10c",
    grid="@earth_relief_01d_g",
    frame=["wSne+tremote dataset", "af"],
)

fig.show()
# fig.savefig(fname="tilemap_lon0to360.png")

Output figure:

[weiji14]

Hmm, a few things:

1. We should probably document that the lonlat range is from -180 to +180, or have a way of converting 0-360 coordinates to -180 to +180 before calling contextily.bounds2img. Specifically these lines
   

   pygmt/pygmt/datasets/tile_map.py

   Lines 27 to 30 in 91e29f1

   	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``.

   
   

   pygmt/pygmt/src/tilemap.py

   Lines 57 to 60 in 91e29f1

   	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``.

2. Ideally, the cross dateline support can be added upstream into contextily (or xyzservices?), and we should at least try to open up an issue to see if they're interested in handling it. I actually need to check if EPSG:3857 and EPSG:4326 tile servers actually support returning cross-dateline tiles. If not, then we may need to add the logic of fetching tiles to the West and East of the international dateline separately and stitch them back together.

[seisman]

2. Ideally, the cross dateline support can be added upstream into contextily (or xyzservices?), and we should at least try to open up an issue to see if they're interested in handling it.

Sounds good.

[seisman]

2. Ideally, the cross dateline support can be added upstream into contextily (or xyzservices?), and we should at least try to open up an issue to see if they're interested in handling it. I actually need to check if EPSG:3857 and EPSG:4326 tile servers actually support returning cross-dateline tiles.

contextily/xyzservices assumes that all tiles come in web mercator (EPSG:3857) (source: geopandas/xyzservices#99 (comment)) and web mercator's center location is (0, 0) and has a bound range (-20037508.34, 20037508.34, -20048966.1, 20048966.1) which corresponds to (-180.0, 180.0, -85.06, 85.06) in WGS84 CRS (source: https://epsg.io/3857). So, I think it's unlikely that xyzservices can return cross-dateline tiles.

If not, then we may need to add the logic of fetching tiles to the West and East of the international dateline separately and stitch them back together.

Yes, we can do it. For example, for a longitude range with west=150, east=200, we can split into two parts, west=150, east=180, and west=-180, east=-160, and call contextily.bounds2img twice, like:

image1, extent1 = contextily.bounds2img(w=150, s=south, e=180, n=north, **kwargs)
image2, extent2 = contextily.bounds2img(w=-180, s=south, e=-160, n=north, **kwargs)

However, the problem is, currently load_tile_map always returns images in EPSG:3857 CRS. We can merge image1 and image2 into a single data array, but we can't merge extent1 and extent2, because extent1 is [16697923.62, 20037508.34, ..., ...] and extent2 is [-20037508.34, -17811118.52, ..., ...].

Since the cross-dateline issue only exists when lonlat=True, I think the solution is,

1. Add a new parameter (e.g., crs) to load_tile_map to allow changing CRS of the returned image (load_map_tiles: Allow a new parameter to control the projection of the returned DataArray #3484)
2. Deal with cross-dateline issue only when lonlat=True and crs is set to something like OGC:CRS84.