From 699befdee093ab4f59bdff96dc48f04eaf55710c Mon Sep 17 00:00:00 2001 From: mkrupczak3 Date: Wed, 6 Nov 2024 11:21:15 -0700 Subject: [PATCH] add DTED Interpolation, remove ray angle rotation (#178) * unfinished DTED interpolation implementation * add IDW for DTED (fix #177), remove earth curvature ray rotation --- .../main/java/com/openathena/DEMParser.java | 85 +++++++++++++++++-- .../main/java/com/openathena/DTEDParser.java | 46 ---------- .../java/com/openathena/MainActivity.java | 5 ++ .../java/com/openathena/TargetGetter.java | 29 +++---- 4 files changed, 98 insertions(+), 67 deletions(-) delete mode 100644 app/src/main/java/com/openathena/DTEDParser.java diff --git a/app/src/main/java/com/openathena/DEMParser.java b/app/src/main/java/com/openathena/DEMParser.java index d066071e..de61364f 100644 --- a/app/src/main/java/com/openathena/DEMParser.java +++ b/app/src/main/java/com/openathena/DEMParser.java @@ -64,6 +64,9 @@ public class DEMParser implements Serializable { try { loadDEM(geofile); } catch (IllegalArgumentException | TiffException e) { + Log.d(TAG, "Could not parse as GeoTIFF file: " + e.getMessage()); + Log.d(TAG, "now trying to parse as DTED file..."); + //e.printStackTrace(); // If GeoTIFF parsing fails, try to parse as DTED try { this.dted = new FileBasedDTED(geofile); @@ -155,8 +158,8 @@ public void loadDEM(File geofile) throws IllegalArgumentException { FileDirectoryEntry fde = directory.get(FieldTagType.GeoKeyDirectory); // Because my sources don't properly store vertical datum, // we have to assume its EGM96 (for GeoTIFFs) and hope it's correct :( - // When we implement the DTED2 and DTED3 format, - // we know it will also be in EGM96 + // For the DTED2 and DTED3 format, + // we know it will be in EGM96 because of the DTED spec verticalDatum = verticalDatumTypes.EGM96; if (fde != null) { @@ -502,16 +505,86 @@ public double getAltFromLatLon(double lat, double lon) throws RequestedValueOOBE } if (this.isDTED) { - Point point = new Point(lat, lon); + // on DEM edge check, if so just return nearest elevation + if (lat == getMaxLat() || lat == getMinLat() || lon == getMaxLon() || lon == getMinLon()) { + Point point = new Point(lat, lon); + try { + double EGM96_altitude = this.dted.getElevation(point).getElevation(); + // DTED vertical datum is height above EGM96 geoid, we must convert it to height above WGS84 ellipsoid + double WGS84_altitude = EGM96_altitude - offsetProvider.getEGM96OffsetAtLatLon(lat,lon); + return WGS84_altitude; + } catch (CorruptTerrainException e) { + throw new CorruptTerrainException("The terrain data in the DTED file is corrupt.", e); + } catch (InvalidValueException e) { + throw new RequestedValueOOBException("getAltFromLatLon arguments out of bounds!", lat, lon); + } + } + + // Determine DTED grid resolution based on level + // For example, Level 0: 1 arc-second, Level 1: 3 arc-seconds, etc. + double gridLatStep = this.dted.getLatitudeInterval(); + double gridLonStep = this.dted.getLongitudeInterval(); + + // Calculate the surrounding grid points + double latFloor = Math.floor(lat / gridLatStep) * gridLatStep; + double lonFloor = Math.floor(lon / gridLonStep) * gridLonStep; + + // Define the four surrounding points + Point p1 = new Point(latFloor, lonFloor); + Point p2 = new Point(latFloor, lonFloor + gridLonStep); + Point p3 = new Point(latFloor + gridLatStep, lonFloor); + Point p4 = new Point(latFloor + gridLatStep, lonFloor + gridLonStep); + try { - double EGM96_altitude = this.dted.getElevation(point).getElevation(); - // DTED vertical datum is height above EGM96 geoid, we must convert it to height above WGS84 ellipsoid - double WGS84_altitude = EGM96_altitude - offsetProvider.getEGM96OffsetAtLatLon(lat,lon); + // Retrieve elevations for the four surrounding points + double e1 = this.dted.getElevation(p1).getElevation(); + double e2 = this.dted.getElevation(p2).getElevation(); + double e3 = this.dted.getElevation(p3).getElevation(); + double e4 = this.dted.getElevation(p4).getElevation(); + + // Create Location objects for each point + Location L1 = new Location(p1.getLatitude(), p1.getLongitude(), e1); + Location L2 = new Location(p2.getLatitude(), p2.getLongitude(), e2); + Location L3 = new Location(p3.getLatitude(), p3.getLongitude(), e3); + Location L4 = new Location(p4.getLatitude(), p4.getLongitude(), e4); + +// Log.d(TAG, "Interpolating DTED DEM using IDW"); +// Log.d(TAG, "Lat Lon: " + lat + "," + lon); +// Log.d(TAG, "L1: " + p1.getLatitude() + "," + p1.getLongitude()); +// Log.d(TAG, "L2: " + p2.getLatitude() + "," + p2.getLongitude()); +// Log.d(TAG, "L3: " + p3.getLatitude() + "," + p3.getLongitude()); +// Log.d(TAG, "L4: " + p4.getLatitude() + "," + p4.getLongitude()); + + // Define the target location + Location target = new Location(lat, lon); + + // Array of neighboring points + Location[] neighbors = new Location[]{L1, L2, L3, L4}; + + // Define the power parameter for IDW + double power = 1.875d; + + // Perform IDW interpolation + double interpolatedAltitude = idwInterpolation(target, neighbors, power); + +// double rawAltitude = this.dted.getElevation(new Point(lat, lon)).getElevation(); +// Log.d(TAG, "rawAltitude: " + rawAltitude); +// Log.d(TAG, "interpolatedAltitude: " + interpolatedAltitude); + + // Convert from EGM96 AMSL orthometric height to WGS84 HAE + double WGS84_altitude = interpolatedAltitude - offsetProvider.getEGM96OffsetAtLatLon(lat, lon); + +// // Convert from EGM96 AMSL orthometric height to WGS84 HAE +// double WGS84_altitude = rawAltitude - offsetProvider.getEGM96OffsetAtLatLon(lat, lon); + + return WGS84_altitude; } catch (CorruptTerrainException e) { throw new CorruptTerrainException("The terrain data in the DTED file is corrupt.", e); } catch (InvalidValueException e) { throw new RequestedValueOOBException("getAltFromLatLon arguments out of bounds!", lat, lon); + } catch (Exception e) { + throw new CorruptTerrainException("An unexpected error occurred while processing DTED data.", e); } } diff --git a/app/src/main/java/com/openathena/DTEDParser.java b/app/src/main/java/com/openathena/DTEDParser.java deleted file mode 100644 index d648a87b..00000000 --- a/app/src/main/java/com/openathena/DTEDParser.java +++ /dev/null @@ -1,46 +0,0 @@ -package com.openathena; - -import com.agilesrc.dem4j.Point; -import com.agilesrc.dem4j.dted.impl.FileBasedDTED; -import com.agilesrc.dem4j.exceptions.CorruptTerrainException; -import com.agilesrc.dem4j.exceptions.InvalidValueException; - -import java.io.File; -import java.io.FileNotFoundException; - -public class DTEDParser { - private FileBasedDTED dted; - private String filePath; - - public DTEDParser(String filePath) throws FileNotFoundException { - this.filePath = filePath; - File file = new File(filePath); - if (!file.exists()) { - throw new FileNotFoundException("The file " + filePath + " does not exist."); - } - try { - this.dted = new FileBasedDTED(file); - } catch (Exception e) { - throw new IllegalArgumentException("Failed to parse the DTED file: " + filePath, e); - } - } - - public double getElevation(double lat, double lon) throws Exception { - Point point = new Point(lat, lon); - try { - return this.dted.getElevation(point).getElevation(); - } catch (CorruptTerrainException e) { - throw new Exception("The terrain data in the DTED file is corrupt.", e); - } catch (InvalidValueException e) { - throw new Exception("The provided latitude and longitude values are invalid.", e); - } - } - - public String getDTEDLevel() { - try { - return this.dted.getDTEDLevel().name(); - } catch (CorruptTerrainException cte) { - return "CORRUPT"; - } - } -} diff --git a/app/src/main/java/com/openathena/MainActivity.java b/app/src/main/java/com/openathena/MainActivity.java index 87a45720..32452721 100644 --- a/app/src/main/java/com/openathena/MainActivity.java +++ b/app/src/main/java/com/openathena/MainActivity.java @@ -851,7 +851,12 @@ public void calculateImage(View view, boolean shouldISendCoT) if (theTGetter != null) { try { + long startTime = System.nanoTime(); result = theTGetter.resolveTarget(y, x, z, azimuth, theta); + long endTime = System.nanoTime(); + long duration = endTime - startTime; + System.out.println("resolveTarget duration in milliseconds: " + duration / 1_000_000.0d); + distance = result[0]; latitude = result[1]; longitude = result[2]; diff --git a/app/src/main/java/com/openathena/TargetGetter.java b/app/src/main/java/com/openathena/TargetGetter.java index 4f260d38..2bda7ed1 100644 --- a/app/src/main/java/com/openathena/TargetGetter.java +++ b/app/src/main/java/com/openathena/TargetGetter.java @@ -116,11 +116,10 @@ public double[] resolveTarget(double lat, double lon, double alt, double azimuth double curAlt = alt; double groundAlt; - // account for curvature of Earth at each Epoch long iterCount = 0; - long iterPerEpoch = 128; - double lastEpochLat = curLat; - double lastEpochLon = curLon; +// long iterPerEpoch = 128; +// double lastEpochLat = curLat; +// double lastEpochLon = curLon; try { groundAlt = myDEMParser.getAltFromLatLon(curLat, curLon); } catch (RequestedValueOOBException e) { @@ -143,17 +142,17 @@ public double[] resolveTarget(double lat, double lon, double alt, double azimuth } altDiff = curAlt - groundAlt; - // account for curvature of Earth at each Epoch - if (iterCount > 0 && iterCount % iterPerEpoch == 0) { - double minArcAngle = minArcAngle(lastEpochLat, lastEpochLon, curLat, curLon); - Log.d(TAG, "iterCount: " + iterCount + ", minArcAngle: " + minArcAngle); - radTheta += Math.toRadians(minArcAngle); // rotate ray downwards by same amount Earth's curvature rotated it upwards - deltaZ = -1.0d * Math.sin(radTheta); - horizScalar = Math.cos(radTheta); - - lastEpochLat = curLat; - lastEpochLon = curLon; - } +// // account for curvature of Earth at each Epoch +// if (iterCount > 0 && iterCount % iterPerEpoch == 0) { +// double minArcAngle = minArcAngle(lastEpochLat, lastEpochLon, curLat, curLon); +// Log.d(TAG, "iterCount: " + iterCount + ", minArcAngle: " + minArcAngle); +// radTheta += Math.toRadians(minArcAngle); // rotate ray downwards by same amount Earth's curvature rotated it upwards +// deltaZ = -1.0d * Math.sin(radTheta); +// horizScalar = Math.cos(radTheta); +// +// lastEpochLat = curLat; +// lastEpochLon = curLon; +// } double avgAlt = curAlt; curAlt += deltaZ;