I want to create a grid of WGS84 ellipsoid (with the distance to the center of the earth).
I found that I could solve this equation where the “radii” is a function of the latitude.
a: mayor axis
b: minor axis
Is there any module that have this formula already incorporate? mapproject?
Or should I do it from scratch?
Seems like a grdmath macro you create and stick in ~/.gmt/grdmath,macros, e.g.
# Accept a b and evaluate ellipsoid on the grid
ELLIPSOID = STO@a POP STO@b POP RCL@a 2 POW Y COSD MUL 2 POW RCL@b 2 POW Y SIND MUL 2 POW ADD etc etc
then
gmt grdmath -Rg -I1 6378137.0 6356752.3142 ELLIPSOID = my.grd
I have many macros for grdmath and gmtmath - perhaps we should show people how to do these things. I am sure only 0.00044% of GMT users are even aware of this capability.
Sure, because compare that with
julia> a = 6378137.0;
julia> b = 6356752.3142;
julia> phi = -90:1:90;
julia> f(phi) = sqrt.( ((a^2 * cosd.(phi)).^2 .+ (b^2 * sind.(phi)).^2) ./ ((a * cosd.(phi)).^2 .+ (b * sind.(phi)).^2) )
f (generic function with 1 method)
julia> col = f(phi);
julia> mat = repeat(col,1,361);
julia> G = mat2grid(mat, hdr=[0,360,-90,90, b, a, 0, 1, 1]);
imshow(G, coast=true)
Never knew this; please do.
Yes, but now we are talking apples and oranges. We can do the same in matlab and python. And the amount of coding in those languages are not necessarily less - usually RPN expressions are more compact.
“More compact” does not necessarily mean “easier to implement” or “more readable/understandable” for each particular user. Often means “more error-prone” and “very steep learning curve”.
I have many macros for grdmath and gmtmath - perhaps we should show people how to do these things.
please do, would be much appreciated.
@mkononets said most of what wanted. The problem with grdmath for equations is it quickly becomes very complex to incomprehensible. I don’t imagine the time it would take me to write that expression.
Fair enough. I speak fluent RPN 
Many thanks to all the answers.
I just wonder if Julia can translate any formula into RPN (if it is good as they said).
Julia is a programming language and like any programming languages it does nothing by itself. It’s people (and now other programs too) who write programs to achieve the wished objectives.
But what would be the interest in translating a formula to RPN when the idea is to execute that formula?
Just curious. But it could be helpfull as an interactive tool to learn RPN.
My impression that grdmath has been created when there was nothing like Python/Julya around with rich native libraries greatly simplifying processing of netcdf grid files and the like.
I am more curious to see what people like Pål do with grids using grdmath rather than exactly how to do it with grdmath.
It was easier than I expected.
ELLIPSOID = STO@a POP STO@b POP RCL@a 2 POW Y COSD MUL 2 POW RCL@b 2 POW Y SIND MUL 2 POW ADD RCL@a Y COSD MUL 2 POW RCL@b Y SIND MUL 2 POW ADD DIV SQRT
then
gmt grdmath -Rg -I1 6378137.0 6356752.3142 ELLIPSOID = my.grd
Now I wonder if any way to extract the values of the elipsoid from GMT. I
I can list the values of the elipsoids with gmt mapproject -Qe. Is any way to take those values as inputs (instead of writing them in the formula)?
ID Date a 1/f
->WGS-84 1984 6378137.000 298.25722356Well, you could have your macro expect a and 1/f and start by computing b = a (1-f) and keep the rest. e.g.
STO@a POP INV 1 SUB NEG RCL@a MUL STO@b POP …
Sorry, my doubt is how to extract the info (values a and 1/f) from the table.
The actual ouptut of gmt mapproject -Qe is:
GMT supports 76 ellipsoids, given below (-> indicates default setting)
ID Date a 1/f
-----------------------------------------------------------
Airy 1830 6377563.396 299.324964600
Airy-Ireland 1830 6377340.189 299.324964600
Andrae 1876 6377104.430 300.000000000
APL4.9 1965 6378137.000 298.250000000
ATS77 1977 6378135.000 298.257000000
Australian 1965 6378160.000 298.250000000
Bessel 1841 6377397.155 299.152812800
Bessel-Namibia 1841 6377483.865 299.152812800
Bessel-NGO1948 1841 6377492.018 299.152810000
Bessel-Schwazeck 1841 6377483.865 299.152812800
Clarke-1858 1858 6378293.639 294.260680000
Clarke-1866 1866 6378206.400 294.978698200
Clarke-1866-Michigan 1866 6378450.047 294.978698200
Clarke-1880 1880 6378249.145 293.465000000
Clarke-1880-Arc1950 1880 6378249.145 293.466307600
Clarke-1880-IGN 1880 6378249.200 293.466021300
Clarke-1880-Jamaica 1880 6378249.136 293.466310000
Clarke-1880-Merchich 1880 6378249.200 293.465980000
Clarke-1880-Palestine 1880 6378300.790 293.466230000
CPM 1799 6375738.700 334.290000000
Delambre 1810 6376428.000 311.500000000
Engelis 1985 6378136.050 298.256600000
Everest-1830 1830 6377276.345 300.801700000
Everest-1830-Kalianpur 1830 6377301.243 300.801740000
Everest-1830-Kertau 1830 6377304.063 300.801700000
Everest-1830-Pakistan 1830 6377309.613 300.801700000
Everest-1830-Timbalai 1830 6377298.556 300.801700000
Fischer-1960 1960 6378166.000 298.300000000
Fischer-1960-SouthAsia 1960 6378155.000 298.300000000
Fischer-1968 1968 6378150.000 298.300000000
FlatEarth 1984 6371008.771 inf
GRS-67 1967 6378160.000 298.247167427
GRS-80 1980 6378137.000 298.257222101
Hayford-1909 1909 6378388.000 297.000000000
Helmert-1906 1906 6378200.000 298.300000000
Hough 1960 6378270.000 297.000000000
Hughes-1980 1980 6378273.000 298.279400000
IAG-75 1975 6378140.000 298.257222000
Indonesian 1974 6378160.000 298.247000000
International-1924 1924 6378388.000 297.000000000
International-1967 1967 6378157.500 298.250000000
Kaula 1961 6378163.000 298.240000000
Krassovsky 1940 6378245.000 298.300000000
Lerch 1979 6378139.000 298.257000000
Maupertius 1738 6397300.000 191.000000000
Mercury-1960 1960 6378166.000 298.300000000
MERIT-83 1983 6378137.000 298.257000000
Modified-Airy 1830 6377340.189 299.324964600
Modified-Fischer-1960 1960 6378155.000 298.300000000
Modified-Mercury-1968 1968 6378150.000 298.300000000
NWL-10D 1972 6378135.000 298.260000000
NWL-9D 1966 6378145.000 298.250000000
OSU86F 1986 6378136.200 298.257220000
OSU91A 1991 6378136.300 298.257220000
Plessis 1817 6376523.000 308.640000000
SGS-85 1985 6378136.000 298.257000000
South-American 1969 6378160.000 298.250000000
Sphere 1984 6371008.771 inf
Struve 1860 6378297.000 294.730000000
TOPEX 1990 6378136.300 298.257000000
Walbeck 1819 6376896.000 302.780000000
War-Office 1926 6378300.583 296.000000000
WGS-60 1960 6378165.000 298.300000000
WGS-66 1966 6378145.000 298.250000000
WGS-72 1972 6378135.000 298.260000000
->WGS-84 1984 6378137.000 298.257223563
Web-Mercator 1984 6378137.000 inf
Moon 2000 1737400.000 inf
Mercury 2000 2439700.000 inf
Venus 2000 6051800.000 inf
Mars 2000 3396190.000 169.894447224
Jupiter 2000 71492000.000 15.414402760
Saturn 2000 60268000.000 10.207994580
Uranus 2000 25559000.000 43.616040956
Neptune 2000 24764000.000 58.543735225
Pluto 2000 1195000.000 infWell, presumably you know the name of the ellipsoid you want, e…, TOPEX so you will have to do
gmt mapproject -Qe 2>&1 | grep TOPEX
since mapproject writes that to stderr…
If you really need to automate ellipsoid building then I guess you could temporarily create a macro in the current directory? e.g.
gmt mapproject -Qe 2>&1 | grep TOPEX | awl '{TOPEX = %lf STO@a POP %lf INV 1 SUB NEG RCL@a MUL …\n", $3, $4} > grdmath.macros
and then just grdmath -R -I TOPEX = topex.grd
Or something like that.
I’m afraid your solution is not correct. Your Earth is thinner at the equator and swelled at the poles when it should be the contrary.
Yes, thanks Joaquim. I have to flip the values of the axis.