I am having issues with the sph2grd module in GMT (6.4.0). I would like to be able to use the spherical harmonic coefficient files available at the ICGEM site (gfc files), but cannot get the sph2grd to generate anything meaningful.
I tested the example in the documentation and the result was the same:
gmt sph2grd @EGM96_to_36.txt -GEGM96_to_36.nc -Rg -I1 -V (see attached)
The exact same result happens with the ICGEM data, i.e. an identical image!
Any ideas as to what is going on?
Also the documentation quotes the use of the -D option for computing the gravity (g) or geoid (n) but this gives an error stating it is not implemented.
The EGM example is very long wavelength so mostly latitudinal changes as you see. Also see cookbook example 39 which runs fine for us and illustrates it is working OK. I cannot speak about gfc files and if they are in the expected format or not though.
I think that the “issue” is that your result is over the sphere, when you should calculate it over the elipsiod in order to get something meaningful. Is it possible?
I stripped off the header parts and removed the column with gfc in it and it generates the same plot as the EGM example. Tried with a file up to D/O 2190 and the same result.
The gfc files are SHM format with an extra column and header section. Example 39 works fine for me also. I have attached a zipped copy of the file in its original unedited format.
This has come out already and my suggestion back then as well as today is: try Geographiclib, sph2grd has several issues. Also, from memory, the Potsdam sites had online tools to generate grids from the spherical harmonics coefficients.
the spherical harmonic coefficient file usually contains the entire field, so the magnitude is relatively large and symmetrical about the equator. You should subtract the normal field before using sph2grd .
You can use GrafLab by Bucha and Janák (2013). It can handle almost all the potential field signals. The results are in XYZ format, which you can use GMT’s surface or nearneighbor to convert into a grid file. Or better still, use the reshape function in MATLAB to convert into a matrix (grid).
Bucha B, Janák J (2013) A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders. Comput Geosci 56:186–196. https://doi.org/10. 1016/j.cageo.2013.03.012