Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

load_map_tiles: Allow a new parameter to control the projection of the returned DataArray #3484

Open
seisman opened this issue Oct 4, 2024 · 0 comments · May be fixed by #3554
Open

load_map_tiles: Allow a new parameter to control the projection of the returned DataArray #3484

seisman opened this issue Oct 4, 2024 · 0 comments · May be fixed by #3554
Labels
enhancement Improving an existing feature
Milestone

Comments

@seisman
Copy link
Member

seisman commented Oct 4, 2024

Originally proposed in #2125 (comment)

Perhaps we should convert the returned image to longitude/latitude because it's more commonly used than the Spherical Mercator coordinate system, following https://contextily.readthedocs.io/en/latest/warping_guide.html.

Most tiles are originally served as Web Mercator (EPSG:3857), and reprojecting to latlon (EPSG:4326) would result in distortion. I'd prefer the reprojection to be a user controlled step (e.g. using rioxarray's .rio.to_crs), because if latlon isn't what the user wants, there would be double reprojection (EPSG:3857 ->EPSG:4326 ->EPSG:????) which is less accurate than single reprojection (EPSG:3857 ->EPSG:????).

I'd prefer the reprojection to be a user controlled step (e.g. using rioxarray's .rio.to_crs)

What about having an option so that users can decide on the desired projection and don't have to learn the syntax of the rioxarray package?

I think we can add a new parameter like crs to control the projection. crs should be default to EPSG:3857. If a different CRS is given, then we can use contextily.warp_tiles to reproject the tiles.

With this new parameter, the following lines in Figure.tilemap can be simplified:

pygmt/pygmt/src/tilemap.py

Lines 141 to 145 in 68a17a0

# Reproject raster from Spherical Mercator (EPSG:3857) to lonlat (OGC:CRS84) if
# bounding box region was provided in lonlat
if lonlat and raster.rio.crs == "EPSG:3857":
raster = raster.rio.reproject(dst_crs="OGC:CRS84")
raster.gmt.gtype = 1 # set to geographic type

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Improving an existing feature
Projects
None yet
1 participant