How to project Cartesian data using PyGMT?

Q&A PyGMT Q&A
1

Greetings dear PyGMT users,

The question is: how can I import a grid in PyGMT which initially is in Mercator CRS and set in the code what is the CRS of the imported data?

It will allow me to plot this grid using frame with geographical coordinates and add other geographic data on top of it.

Thank you.

2

Hi @Igor,

Are you working with NetCDF grids? GeoTiffs? Or another format? Is it a single band or multi-band grid?

If you want to load it into an xarray.DataArray, what I would recommend is rioxarray's reproject function (see https://corteva.github.io/rioxarray/stable/examples/reproject.html) or rasterio’s reproject to do the CRS conversion (both based on GDAL). You can then use PyGMT to plot the reprojected xarray grid later (and overlay other map layers on top).

If you’re after a pure GMT way, then grdproject is what you’re after.

1 Like
3

@weiji14, thank you very much for your answer and a list of possible solutions.

I am working with both surfer grids and xyz data.

Due to inevitable loss of grid information while reprojecting and resampling, I decided finally to keep my grid in Mercator projection, and project other data on top of it. These tutorials have helped me out:

  1. https://docs.generic-mapping-tools.org/dev/gallery/ex27.html#example-27
  2. https://docs.generic-mapping-tools.org/dev/gallery/ex28.html#example-28

Thanks a lot for your reply once again.

1 Like