Plotting spherical coordinates with logarithmic scale

Added by Lars O 8 months ago

Dear all,

I just started again to (ab)use GMT as an interpolation and plotting tool for large, irregular data-sets (approx. 150000 data-points, too much for most other tools I have at hand). This is not geographic data, but similar enough in structure to get it gridded and plotted. However, there is some fine-tuning needed and I am somehow not getting any further. My data is in three columns θ,φ,val (with φ being the azimuth). The val-range is typically over 5-8 decades, so I want to plot it with a logarithmic scale. Center of my plots will always be θ=φ=0. Here is what I do so far:

1) generate a triangle network and a grid, swapping the first two columns (since GMT expects the azimuth first):

gmt triangulate test.in -R-180/180/0/90 -Gtest.grd -I.1 -M -Z -:i > test.out

2) build a color-map from the mesh (I also tried makecpt, but this seams easier for now):

gmt grd2cpt test.grd -Cjet -Qo -Z > jetlog.cpt

3) build a shading layer from the mesh:

gmt grdgradiant test.grd -Gtest.shd -A90/50 -Ne.6 -fg

4) make a polar cylindrical plot:

gmt grdimage test.grd -JP10c -Cjetlog.cpt -Bg10a0p/g10a10p -Itest.shd > test2d.ps

5) make a perspective view of the plot

gmt grdview test.grd -JP10c -Cjetlog.cpt -Itest.shd > test3d.ps

This leads me to two questions / problems:

The result of stept 4) is quite close to what I want. The colorscale is however not distributed over the full z-range - I get a more or less uniform blue circle with only the tiny peak regions get the color-range. This is not really what I expected from a log colormap. Generating it with -Qi (in step 2) gives a niver output, but seams to be illogical since my input grid is not log-scaled.

Step 5) gives me only a flat wireframe mesh right now. Is it actually possible to generate a relief like a height-map of the polar projection? Or is this just beyond what this tool is ment for?

I went through a bunch of manpages and tutorials now, but could not find an obvious answer probably since my application is a bit beyond the typical use of GMT. However it would be great if someone can hint me at a possible solution, since GMT is close to become a very valuable plotting tool for me.

Regards,
Lars.


Replies (4)

RE: Plotting spherical coordinates with logarithmic scale - Added by Paul 8 months ago

Suspect the problem lies in step 2 which is doing a histogram-equalized CPT where all nodes are given equal weight (i.e., a Cartesian grid weighting). However, that is not true in your polar projection where lots of points represent much smaller areas. I would try to run grdproject on your grid using your map projection, which will create a Cartesian grid in the polar projection. Now try grd2cpt on that grid to get another cpt and see if this helps.

RE: Plotting spherical coordinates with logarithmic scale - Added by Lars O 5 months ago

Hi Paul,

thanks to your hints and the huge amount documentation and supporting material available, I have managed to get the expected plot. I run a sequence of triangulate, makecpt, grdgradient and grdimage to get the polar projection of my hemishperical distribution of scattered data-points:

gmt makecpt -Cjet -T${SCALEMAX}/${SCALEMIN}/1 -Qi -Z > ${WDIR}/test.cpt
gmt triangulate test.in -R-180/180/0/90 -Gtest.grd -I.1 -M -Z -:i > test.out
gmt grdgradient test.grd -Gtest.grad=nb/a -E30/45/0.55/0.6/0.4/10 -Ne0.6 -fg
gmt grdimage test.grd -JP10c/-90 -Ctest.cpt -Bg10a30p/g10a10p -Itest.grad > test.ps

SCALEMAX and SCALEMIN are pre-calculated limits. This gives me a result su as the attached plot:

The one problem I still have is that there is a gap between the azimuth -180° and +180° leading to a grey line (downward in the illustration). Is there any possible solution a close that gap in my plots?

Best regards, Lars.

model_t_020_0900.jpg - Polar projection of scattered hemispherical data-points. (196 KB)

model_t_020_0900.jpg (196 KB)

RE: Plotting spherical coordinates with logarithmic scale - Added by Lars O 2 months ago

Hi all,

I just wanted to share what I found out recently. The problem is that GMT (5) seams to be able to handle periodicity only if the region of interest range starts at 0. Therefore, my longitutes from -180 to 180 were not connected. Rewriting the negative longitudes so that the range was all positive 0 to 360 gave me the continous output I had been aiming at.

Cheers, Lars.

RE: Plotting spherical coordinates with logarithmic scale - Added by Lars O 2 months ago

Hi all,

I am following up to this old thread rather than opening a new one, since the application is exotic enough and I would like to keep this together.

Now that I have solved 2d plotting, I tried to generate a elevation grid with grdview.

gmt grdview test.grd -JP10c/-90 -Bg10a30p/g10a10p -Ctest.cpt -Jz.5 -pz130/30 -P -Rtest.grd -f0x,1y -nb > test.ps

However, I am facing two new issues now:

  • The interpolation connects some vertices where it shouldn't. Basically I would like to keep interpolation bound to x,y.
  • The periodical data, which I got "connected" at azimuth 0°, 360°, is ignored by grdview - there is a discontinuity at azimuth 0°.
    !!
    I am still working with GMT 5.x - should I upgrade to 6.x, or are these problems not solved yet?

Best regards, Lars.
Result of grdview for polar projection.

test.jpg - Result of grdview for polar projection. (207 KB)

(1-4/4)