[PROJ] Read vertical deflection from egm (96, 08, etc.)

[Charles Karney]

There are two problems here:

(1) you are numerically differentiating a noisy function, the height of
the geoid found by interpolating into a grid of values.

(2) the geoid is a surface of constant gravitional potential, however
the gravity vector is not normal to it!  This is because the geoid goes
under the land mass of the earth.

I recommend using the full gravity model to determine the vertical
deflection.  GeographicLib's GravityModel class gives you the tools to
do this.

Note: GeographicLib used to allow you to compute the gradient of the
geoid along the lines you suggest.  However, I took this functionality
out in 2016 (version 1.46) because of the shortcomings listed above.

   --Charles

On 6/26/19 10:13 AM, Fabian Gross wrote:
> Dear proj-members,
> 
> is there a way to read the vertical deflection (xi, eta) directly with 
> proj?
> 
> My work-around is sampling the geoid-grid shift:
> 
> ```python
> *def *get_vertical_deflection(lon, lat):
> /"""
>     take origin and 1 east + 1 north on GEOID, transform to topocentric 
> and check
> /*:param*/lon: longitude of point of interest on WGS84 (°)
> /*:param*/lat: latitude of point of interest on WGS84 (°)
> /*:return*/: (xi, eta) deflection in meridional and prime vertical (rad)
>     """
> /eps = 0.0001 /# °
> /lons = np.array([lon, lon + eps, lon])
>     lats = np.array([lat, lat,       lat + eps])
>     geoid_heights = get_geoid_height_egm96(lons, lats)
> 
> /# W. E. Featherstone, "The Use and Abuse of Vertical Deflection" p. 4
> /e = eccentricity(f_wgs84)
>     phi = np.deg2rad(lat)
> 
> /# meridional, north/south
>     # R_M, should be radius of curvature in the meridian at point of 
> interest
> /rad_meridian = get_radius_of_curvature_meridian(a_wgs84, e, phi)
> /# d_North = R d_lat; in rad
> /xi = -(geoid_heights[2] - geoid_heights[0])/(rad_meridian * 
> np.deg2rad(eps))
> 
> /# prime vertical, east/west
>     # R_N, should be radius of curvature in the prime vertical at point 
> of interest
> /rad_prime_vertical = get_radius_of_curvature_prime_vertical(a_wgs84, e, 
> phi)
> /# d_East = R cos(lat) d_long; in rad
> /eta = -(geoid_heights[1] - geoid_heights[0])/(rad_prime_vertical * 
> np.deg2rad(eps) * np.cos(phi))
> *return *xi, eta
> 
> ```
> Over which (arc) distance should I sample the geoid gtx?
> 
> Kind regards
> 
> Fabian Gross
> Telefon +49 (711) 648 71-995
> _________________________________________________*
> sbp* *
> schlaich
> bergermann partner*
> 
> Beratende Ingenieure
> für erneuerbare Energie
> 
> Stuttgart . Berlin . New York
> São Paulo . Shanghai . Paris
> *
> sbp sonne gmbh*
> 
> Markus Balz Dipl.Ing. (FH)
> Andreas Keil Dipl.Ing.
> 
> Schwabstrasse 43
> 70197 Stuttgart
> Telefon +49 (711) 648 71-0
> _
> __www.sbp.de_ <http://www.sbp.de/>
> _________________________________________________
> 
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
>