Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[mock_uss, monitorlib] Add EGM96 datum conversion #342

Merged
merged 4 commits into from
Nov 15, 2023

Conversation

BenjaminPelletier
Copy link
Member

@BenjaminPelletier BenjaminPelletier commented Nov 14, 2023

U-space requires USSPs to deliver messages with MSL altitude (using EGM96 or EGM2008 geoid as the datum) to authorised users while ASTM F3411-22a transmits altitude between USSPs using the WGS84 ellipsoid as the datum. In preparation to incorporate MSL altitude observations, this PR implements an algorithm to determine the EGM96 geoid offset from the WGS84 ellipsoid and uses it to compute MSL altitude in mock_uss.

The computation is performed as part of the CI (during the NetRID tests). To verify accuracy, I used the following code:

import matplotlib.pyplot as plt
import numpy as np
from s2sphere import LatLng
from monitoring.monitorlib.geo import egm96_geoid_offset

for p in (LatLng.from_degrees(34, -118), ):
    dh = egm96_geoid_offset(p)
    print(f"{p.lat().degrees}, {p.lng().degrees}: {dh:.2f}")

lats = np.arange(-90, 90.1, 0.25)
lngs = np.arange(0, 360, 0.25)
geoid = np.zeros((lats.size, lngs.size), dtype=float)
for i, lat in enumerate(lats):
    for j, lng in enumerate(lngs):
        geoid[i, j] = egm96_geoid_offset(LatLng.from_degrees(-lat, lng))
plt.imshow(geoid, cmap="hot", extent=(0, 360, -90, 90))
plt.xlabel("Longitude (degrees)")
plt.ylabel("Latitude (degrees)")
cbar = plt.colorbar()
cbar.set_label("EGM96 geoid height above WGS84 ellipsoid (meters)", rotation=270)
plt.savefig("egm96.png")

The absolute accuracy of one point was manually verified (34,-118 -> -34.16 meters, versus this calculator), and then the overall shape was verified visually using this output:

egm96

@BenjaminPelletier BenjaminPelletier marked this pull request as ready for review November 14, 2023 19:50
Copy link
Contributor

@barroco barroco left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me 👍 Could you please add the scale on the plot you attached to the PR to help readers evaluate the accuracy?

grid_path = os.path.join(os.path.dirname(__file__), "assets/WW15MGH.DAC")
grid = (
np.fromfile(grid_path, ">i2").reshape(
int(180 / grid_size + 1), int(360 / grid_size)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is not very clear why you added a buffer to the latitude argument and not for longitude, could you please clarify and perhaps add it in comments?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've switched to use enumerated latitudes and longitudes rather than counting them separately here.

/ 100
)
_egm96 = Spline(
np.arange(-90, 90 + grid_size / 2, grid_size),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess the previous response will clarify the latitude buffer here.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The latitude enumeration now has a comment that it's on the [-90, 90] interval. arange does not have an option for an inclusive endpoint, so the targeted endpoint must extend further than the desired endpoint, but not so far as to reach the next step. Adding one half-step works for this.

lat = p.lat().degrees
if lat < -90 or lat > 90:
raise ValueError(f"Cannot compute EGM96 geoid offset at latitude {lat} degrees")
return _egm96.ev(-lat, lng)
Copy link
Contributor

@barroco barroco Nov 15, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you please clarify in a comment why latitude needs to be inverted here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added

@BenjaminPelletier
Copy link
Member Author

Could you please add the scale on the plot you attached to the PR to help readers evaluate the accuracy?

Done; comment updated.

@BenjaminPelletier BenjaminPelletier merged commit 2825146 into interuss:main Nov 15, 2023
8 checks passed
@BenjaminPelletier BenjaminPelletier deleted the egm96 branch November 15, 2023 22:28
github-actions bot added a commit that referenced this pull request Nov 15, 2023
* Add EGM96 datum conversion

* make format

* Address comments and `make format` 2825146
github-actions bot added a commit to openskies-sh/monitoring that referenced this pull request Nov 16, 2023
* Add EGM96 datum conversion

* make format

* Address comments and `make format` 2825146
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants