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

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Fabian Gross
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
_________________________________________________

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: 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
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/proj
>
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj