Problem with sph2grd

Q&A PyGMT Q&A
1

Hello,

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.

Thanks

Lester

Capture_sph2grd

2

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.

3

Any ideas as to what is going on?

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?

4

I have used the the gfc files. The format is ok.

5

Hi Paul,

Here is a brief snippet of the ICGEM gfc file. The data are a global gravity dataset from GOCE

begin_of_head ======================================================================================
modelname IGGT_R1
product_type gravity_field
earth_gravity_constant 0.3986004415E+15
radius 0.6378136460E+07
max_degree 240
norm fully_normalized
tide_system tide_free
errors formal

key L M C S sigma C sigma S
end_of_head ========================================================================================
gfc 0 0 0.100000000000D+01 0.000000000000D+00 0.000000000000D+00 0.000000000000D+00
gfc 1 0 0.000000000000D+00 0.000000000000D+00 0.000000000000D+00 0.000000000000D+00
gfc 2 0 -0.484165361814D-03 0.000000000000D+00 0.967844555286D-10 0.000000000000D+00
gfc 3 0 0.957268237763D-06 0.000000000000D+00 0.615901145325D-10 0.000000000000D+00
gfc 4 0 0.539995010401D-06 0.000000000000D+00 0.569504078188D-10 0.000000000000D+00
gfc 5 0 0.686749531884D-07 0.000000000000D+00 0.420682043944D-10 0.000000000000D+00
gfc 6 0 -0.149943612057D-06 0.000000000000D+00 0.338700680041D-10 0.000000000000D+00
gfc 7 0 0.904990598750D-07 0.000000000000D+00 0.271732517185D-10 0.000000000000D+00
gfc 8 0 0.494953535063D-07 0.000000000000D+00 0.229934978129D-10 0.000000000000D+00
gfc 9 0 0.280091597997D-07 0.000000000000D+00 0.193642902066D-10 0.000000000000D+00

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.

Lester

IGGT_R1.zip (845.8 KB)

SourceICGEM data source

6

Calculating over the Ellipsoid using -E gives same result although it does say it’s not implemented, same as -D

7

Exactly, that is what I suspect. By default is calculated over the sphere and it seems that is the only option.

8

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.

9

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 .

10

How does one remove normal field from the spherical harmonic coefficient file? Not too familiar working directly with this data.

11

You can find it in Heiskanen and Moritz’s Physical Geodesy .
By the way, -D option is not yet implemented in GMT. You can use the calculation service provided by ICGEM. ICGEM International Center for Global Gravity Field Models (gfz-potsdam.de)

12

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