Plotting a projected grid on a map with the same projection

Q&A
1

I have a grid that is in cartesian coordinates (km) calculated from a Lambert equal-area azimuthal projection centred on Iceland (-JA-18/65/x -C-18/65). I would like to plot it on a geographical map using the same projection.

If I create a map e.g. using gmt basemap -JA-18/65/20.5/8, I get a circular map that goes up to 20.5° away from the centre of the projection. A simple way to plot my grid on this would be to use grdproject -I to project my grid into degrees, and then grdimage. This seems suboptimal to me because it just means doing one projection followed by the exact inverse (because grdimage is projecting the geographical grid back into the Lambert grid). Moreover, this results in reduced resolution in the lower latitudes when projecting to degrees. I suppose I could get around this by resampling the original grid at smaller intervals but this doesn’t feel like a very satisfying solution.

I assume another way to do this would be to use grdimage -JX with the original grid, and somehow get the dimensions of the plot to match up — but this comes with the added challenge of how to cut off the grid at the circular map boundary. Or should I just give up on my idea of having a circular map, as it is too complicated? (And I can project everything else I want to plot before plotting and just use -JX the whole time?)

Does anyone have a simple solution to this? Have I missed some option where you can plot a projected grid straight onto a map with the same projection?

Thanks

2

Hello @tmerry1. Welcome to the GMT forum!!

I think that you just need to plot your grid with cartesian coordinates with -JX.
I think that the example 28 could help you.

Federico

3

Thanks Federico! I guess that you are probably right.

I also maintained the circular plot by first using mapproject -Wm to find the radius of my plot in km for the -JA-18/65/20.5/x projection (given by the output ymin and ymax, i.e. third and fourth numbers). So then I could generate a circular polygon with that radius to use as input for gmt clip when plotting my results with -JX -R[ymin]/[ymax]/[ymin]/[ymax] (allowing me to switch between the two projections to plot on the same map).