EGM2008: added undulation data and bilinear computation.

Includes data structure for caching. Part of #224.
2021-03-28 18:53:37 +02:00
parent b395c19832
commit fa50db08c3
3 changed files with 38419 additions and 0 deletions
@@ -0,0 +1,229 @@
+package de.dennisguse.opentracks.util;
+
+import android.content.Context;
+
+import androidx.test.core.app.ApplicationProvider;
+
+import org.junit.Test;
+import org.junit.runner.RunWith;
+import org.junit.runners.JUnit4;
+
+import java.io.DataInputStream;
+import java.io.EOFException;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.InputStreamReader;
+import java.nio.charset.StandardCharsets;
+import java.time.Instant;
+
+import de.dennisguse.opentracks.content.data.TrackPoint;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertNotEquals;
+
+@RunWith(JUnit4.class)
+public class EGM2008UtilsTest {
+
+    private static final double MAX_BILINEAR_ERROR = 0.478;
+
+    private final Context context = ApplicationProvider.getApplicationContext();
+
+    @Test
+    public void fileVerification() throws IOException {
+        // given
+        int expectedLength = 18671444;
+        int expectedHeaderLength = 404;
+
+        String expectedHeader = "P5\n" +
+                "# Geoid file in PGM format for the GeographicLib::Geoid class\n" +
+                "# Description WGS84 EGM2008, 5-minute grid\n" +
+                "# URL http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008\n" +
+                "# DateTime 2009-08-29 18:45:00\n" +
+                "# MaxBilinearError 0.478\n" +
+                "# RMSBilinearError 0.012\n" +
+                "# MaxCubicError 0.294\n" +
+                "# RMSCubicError 0.005\n" +
+                "# Offset -108\n" +
+                "# Scale 0.003\n" +
+                "# Origin 90N 0E\n" +
+                "# AREA_OR_POINT Point\n" +
+                "# Vertical_Datum WGS84\n" +
+                "4320    2161\n" +
+                "65535" +
+                "\n";
+
+        // when
+        try (InputStream inputStream = context.getResources().openRawResource(EGM2008Utils.EGM2008_5_DATA)) {
+            assertEquals(expectedLength, inputStream.available());
+
+            InputStreamReader reader = new InputStreamReader(inputStream, StandardCharsets.US_ASCII);
+            char[] data = new char[expectedHeaderLength];
+            int length = reader.read(data);
+
+            // then
+            assertEquals(expectedHeaderLength, length);
+            assertEquals(expectedHeader, new String(data));
+        }
+    }
+
+    @Test
+    public void data_Northpole() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(90);
+        trackPoint.setLongitude(0.1);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(-14.8980, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void data_Southpole() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(-90);
+        trackPoint.setLongitude(0);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(30.15, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void data_Southpole2() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(-90);
+        trackPoint.setLongitude(180);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(30.15, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void data_Southpole3() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(-90);
+        trackPoint.setLongitude(-180);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(30.15, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void data_0() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(0);
+        trackPoint.setLongitude(0);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(-17.2260, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void data_Berlin() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(52.530644);
+        trackPoint.setLongitude(13.383068);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(-39.4865, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void data_Berlin_Caching() throws IOException {
+        // given
+        TrackPoint trackPoint1 = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint1.setLatitude(52.530644);
+        trackPoint1.setLongitude(13.383068);
+        trackPoint1.setAltitude(0);
+
+        TrackPoint trackPoint2 = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint2.setLatitude(52.530000);
+        trackPoint2.setLongitude(13.380000);
+        trackPoint2.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint1.getLocation());
+
+        // then
+        assertNotEquals(altitude_egm2008.correctAltitude(trackPoint1.getLocation()), altitude_egm2008.correctAltitude(trackPoint2.getLocation()), 0.0001);
+    }
+
+    @Test
+    public void data_MaxUndulation() throws IOException {
+        // given
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(-8.417);
+        trackPoint.setLongitude(147.367);
+        trackPoint.setAltitude(0);
+
+        // when
+        EGM2008Utils.EGM2008Correction altitude_egm2008 = EGM2008Utils.createEGM2008Correction(context, trackPoint.getLocation());
+
+        // then
+        assertEquals(-85.824, altitude_egm2008.correctAltitude(trackPoint.getLocation()), MAX_BILINEAR_ERROR);
+    }
+
+    @Test
+    public void getIndices() {
+        TrackPoint trackPoint = new TrackPoint(TrackPoint.Type.TRACKPOINT, Instant.ofEpochMilli(0));
+        trackPoint.setLatitude(90);
+        trackPoint.setLongitude(0);
+        assertEquals(new EGM2008Utils.Indices(0, 0), EGM2008Utils.getIndices(trackPoint.getLocation()));
+
+        trackPoint.setLatitude(0);
+        trackPoint.setLongitude(0);
+        assertEquals(new EGM2008Utils.Indices(1080, 0), EGM2008Utils.getIndices(trackPoint.getLocation()));
+
+        trackPoint.setLatitude(-90);
+        trackPoint.setLongitude(180);
+        assertEquals(new EGM2008Utils.Indices(2160, 4320 / 2), EGM2008Utils.getIndices(trackPoint.getLocation()));
+
+        trackPoint.setLatitude(-90);
+        trackPoint.setLongitude(-180);
+        assertEquals(new EGM2008Utils.Indices(2160, 0), EGM2008Utils.getIndices(trackPoint.getLocation()));
+    }
+
+    @Test
+    public void getUndulationRaw_ok() throws IOException {
+        try (DataInputStream dataInputStream = new DataInputStream(context.getResources().openRawResource(EGM2008Utils.EGM2008_5_DATA))) {
+            assertEquals(40966, EGM2008Utils.getUndulationRaw(dataInputStream, new EGM2008Utils.Indices(0, 0)));
+            assertEquals(41742, EGM2008Utils.getUndulationRaw(dataInputStream, new EGM2008Utils.Indices(1081, 0)));
+            assertEquals(25950, EGM2008Utils.getUndulationRaw(dataInputStream, new EGM2008Utils.Indices(2160, 4319)));
+        }
+    }
+
+    @Test(expected = EOFException.class)
+    public void getUndulationRaw_error() throws IOException {
+        try (DataInputStream dataInputStream = new DataInputStream(context.getResources().openRawResource(EGM2008Utils.EGM2008_5_DATA))) {
+            assertEquals(0, EGM2008Utils.getUndulationRaw(dataInputStream, new EGM2008Utils.Indices(2161, 4320)), 0.01);
+        }
+    }
+}
@@ -0,0 +1,176 @@
+package de.dennisguse.opentracks.util;
+
+import android.content.Context;
+import android.location.Location;
+
+import androidx.annotation.NonNull;
+import androidx.annotation.VisibleForTesting;
+
+import java.io.DataInputStream;
+import java.io.IOException;
+import java.util.Objects;
+
+import de.dennisguse.opentracks.R;
+
+/**
+ * Converts WGS84 altitude to EGM2008 (should be close to height above sea level).
+ * <p>
+ * Uses <a href="https://geographiclib.sourceforge.io/">GeographicLib</a>] EGM2008 5minute undulation data.
+ * https://geographiclib.sourceforge.io/html/geoid.html
+ * <p>
+ * File starts at 90N, 0E (North pole) and is encoded in parallel bands as unsigned shorts.
+ */
+public class EGM2008Utils {
+
+    static final int EGM2008_5_DATA = R.raw.egm2008_5;
+
+    private static final int HEADER_LENGTH = 404;
+
+    private static final int RESOLUTION_IN_MINUTES = 60 / 5;
+
+    private static final int LATITUDE_CORRECTION = 360 * RESOLUTION_IN_MINUTES;
+
+    private EGM2008Utils() {
+    }
+
+    public static EGM2008Correction createEGM2008Correction(Context context, Location location) throws IOException {
+        Indices indices = getIndices(location);
+
+        try (DataInputStream dataInputStream = new DataInputStream(context.getResources().openRawResource(EGM2008_5_DATA))) {
+            return new EGM2008Correction(indices, dataInputStream);
+        }
+    }
+
+    @VisibleForTesting
+    static int getUndulationRaw(DataInputStream dataInputStream, Indices indices) throws IOException {
+        dataInputStream.reset();
+        int absoluteIndex = indices.getAbsoluteIndex();
+        return getUndulationRaw(dataInputStream, absoluteIndex);
+    }
+
+    private static int getUndulationRaw(DataInputStream dataInputStream, int undulationIndex) throws IOException {
+        dataInputStream.reset();
+        int index = HEADER_LENGTH + undulationIndex * 2;  //byte size is 2
+        long ignored = dataInputStream.skip(index);
+
+        return dataInputStream.readUnsignedShort();
+    }
+
+    @VisibleForTesting
+    static Indices getIndices(Location location) {
+        double latitude = -location.getLatitude() + 90;
+        int latitudeIndex = (int) (latitude * RESOLUTION_IN_MINUTES);
+
+        double longitude;
+        if (location.getLongitude() >= 0) {
+            longitude = location.getLongitude();
+        } else {
+            longitude = 180 + Math.abs(location.getLongitude());
+        }
+        int longitudeIndex = (int) (longitude * RESOLUTION_IN_MINUTES);
+
+        if (longitudeIndex >= 360 * RESOLUTION_IN_MINUTES) {
+            longitudeIndex = 0;
+        }
+        return new Indices(latitudeIndex, longitudeIndex);
+    }
+
+    public static class EGM2008Correction {
+
+        protected final Indices indices;
+
+        protected final int v00;
+        protected final int v10;
+        protected final int v01;
+        protected final int v11;
+
+        public EGM2008Correction(Indices indices, DataInputStream dataInputStream) throws IOException {
+            this.indices = indices;
+
+            v00 = getUndulationRaw(dataInputStream, indices);
+            if (!isSouthPole()) {
+                v10 = getUndulationRaw(dataInputStream, indices.offset(0, 1));
+                v01 = getUndulationRaw(dataInputStream, indices.offset(1, 0));
+                v11 = getUndulationRaw(dataInputStream, indices.offset(1, 1));
+            } else {
+                v10 = 0;
+                v01 = 0;
+                v11 = 0;
+            }
+        }
+
+        public boolean canCorrect(@NonNull Location location) {
+            return indices.getAbsoluteIndex() == getIndices(location).getAbsoluteIndex();
+        }
+
+        public double correctAltitude(@NonNull Location location) {
+            if (!canCorrect(location))
+                throw new RuntimeException("Undulation data not loaded for this location.");
+            if (!location.hasAltitude())
+                throw new RuntimeException("Location has no altitude");
+
+            double undulationRaw;
+            if (isSouthPole()) {
+                // No bilinear interpolation on South Pole (not worth the time)
+                undulationRaw = v00;
+            } else {
+                double fLongitude = location.getLongitude() * RESOLUTION_IN_MINUTES -
+                        (int) (location.getLongitude() * RESOLUTION_IN_MINUTES);
+
+                double fLatitude = (-location.getLatitude() + 90) * RESOLUTION_IN_MINUTES
+                        - (int) ((-location.getLatitude() + 90) * RESOLUTION_IN_MINUTES);
+
+                // Bilinear interpolation (optimized; taken from GeopgrahicLib/Geoid.cpp)
+                double
+                        a = (1 - fLongitude) * v00 + fLongitude * v01,
+                        b = (1 - fLongitude) * v10 + fLongitude * v11;
+                undulationRaw = (1 - fLatitude) * a + fLatitude * b;
+            }
+            // Bilinear interpolation (not optimized)
+            //undulationRaw = v00 * (1 - fLongitude) * (1 - fLatitude)
+            //        + v10 * fLongitude * (1 - fLatitude)
+            //        + v01 * (1 - fLongitude) * fLatitude
+            //        + v11 * fLongitude * fLatitude;
+
+            double h = 0.003 * undulationRaw - 108;
+
+            return location.getAltitude() - h;
+        }
+
+        private boolean isSouthPole() {
+            return indices.latitudeIndex == 2160;
+        }
+    }
+
+    @VisibleForTesting
+    static class Indices {
+        final int latitudeIndex;
+        final int longitudeIndex;
+
+        Indices(int latitudeIndex, int longitudeIndex) {
+            this.latitudeIndex = latitudeIndex;
+            this.longitudeIndex = longitudeIndex;
+        }
+
+        Indices offset(int latitudeOffset, int longitudeOffset) {
+            return new Indices(latitudeIndex + latitudeOffset, longitudeIndex + longitudeOffset);
+        }
+
+        int getAbsoluteIndex() {
+            return latitudeIndex * LATITUDE_CORRECTION + longitudeIndex;
+        }
+
+        @Override
+        public boolean equals(Object o) {
+            if (o == null || getClass() != o.getClass()) return false;
+            Indices indices = (Indices) o;
+            return latitudeIndex == indices.latitudeIndex &&
+                    longitudeIndex == indices.longitudeIndex;
+        }
+
+        @Override
+        public int hashCode() {
+            return Objects.hash(latitudeIndex, longitudeIndex);
+        }
+    }
+}