# [PROJ] Read vertical deflection from egm (96, 08, etc.) Classic List Threaded 2 messages Reply | Threaded
Open this post in threaded view
|

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

 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 - geoid_heights)/(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 - geoid_heights)/(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.)

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