General method for mapping polar region data

Q&A
1

Hi all,

I have some geotiff data for the polar region that are already projected (i.e., in the polar stereographic projection coordinate system, the units of the coordinates XY of a point are meters).
I want to map them with GMT using the -Js (General stereographic) projection.

Since Js/JS only accepts data in geographic coordinates, I need to reproject this data. However, as the actual resolution varies for each latitude, using a fixed lower ppd value for the reprojection reduces the resolution of the data, and using a higher ppd value is time consuming.

I would like to know if there is a general process for handling polar region data. For data in the projected coordinate systems, if one wants to use a projection like Js to map the data, what should one do?

I tried to use grdimage -Jx first, mapping the data under the Cartesian coordinate system, and then use basemap -Js to draw the geographic coordinates at the same scale, but this requires manually aligning the origin of the two coordinate systems. Is there any way to anchor the origin of the coordinate system to a fixed position on the paper?

Thank you!

2

mapproject is what you’re looking for :slight_smile:

here I used grdproject because I had a gridded file. But I guess it can still serve as a support to complete the documentation.

3

Thanks a lot, I tried it and the result looks good. From the parameters it looks like that grdproject works similarly to gdalwarp, still using the same ppd resolution for different latitudes, so I didn’t try before asking here. But specifying a higher resolution produces a netcdf result that will be smaller than the geotiff result produced by gdalwarp. This is basically what I want.

However, I encountered a problem when I used a custom ellipsoid. Using --PROJ_ELLIPSOID=6378137 to specify an ellipsoid would produce the following error:

ERROR 1: PROJ: proj_create: Error 1027 (Invalid value for an argument): pj_init_ctx: Must specify ellipsoid or sphere
grdproject [WARNING]: gmtnc_grd_info failed to convert the proj4 string
+proj=longlat +ellps=unnamed +no_defs
  to WKT

But this problem does not seem to affect the results, the data can still be mapped correctly. This may be a bug, I’m considering whether to open an issue. Anyway, the problem is basically solved, thank you again!