How to grid data for polar plots

Added by Paul 8 months ago

I have a smooth function of (theta, phi) in spherical coordinates and I want to plot it with a polar stereographic plot with GMT. Currently I have a C program which defines a regularly spaced lat/lon grid with 1024 longitude points and 512 latitude points. It then just makes an XYZ file with columns "lon lat value". I convert this XYZ file to a GRD file with this command:

gmt xyz2grd ${XYZFILE} -G${GRDFILE} -Rd -I1024+/512+

and then plot it with grdimage. This generally works ok, however I believe the polar view is not the best it could be, since on a regular lat/lon grid the points will bunch up near the poles. I can see some peculiar features near the north pole, and I want to be sure its not a grid spacing effect and is really in the function I am plotting. I am interested in using a more sophisticated method for generating my grid points on the sphere, such as that described here:

https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf

or the discussion here:

https://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere

or also here:

https://www.sciencedirect.com/science/article/pii/S0010465509001180

These algorithms appear to avoid the problem of having grid points bunched up near the poles, but some of the methods result in a non-regular grid spacing. My question is: can GMT handle such grids? If I generate such a grid in XYZ format, and then use xyz2grd on it, can I plot the result using grdimage as usual? The grdimage program requires knowledge of the grid spacing with the -I option, so how do you get around this for these more sophisticated non-regularly spaced spherical grids?

Does anyone know of an example script somewhere which plots one of these non-regular grids?


Replies (1)

RE: How to grid data for polar plots - Added by Paul 8 months ago

GMT does not accept non-equidistant grids. Since grdimage results in a PostScript image which by definition is equidistant there has to be resampling at some step, even if GMT could ingest a non-regular grid. Alternatively, if you have (x,y,z) triplets you can plot the surface directly with pscontour, in which case there are no interpolations and the surface is represented by Delaunay triangles. Finally, you could generate your points onto an equidistant polar grid directly and plot this as a Cartesian grid. Now there would no resampling, and you could overlay the polar basemap separately. This is similar to how UTM grids are plotted.

(1-1/1)