Skip to content

Commit

Permalink
add DTED Interpolation, remove ray angle rotation (#178)
Browse files Browse the repository at this point in the history 
* unfinished DTED interpolation implementation

* add IDW for DTED (fix #177), remove earth curvature ray rotation
  • Loading branch information
@mkrupczak3
mkrupczak3 authored Nov 6, 2024
1 parent 5058cba commit 699befd
Show file tree
Hide file tree
Showing 4 changed files with 98 additions and 67 deletions.
85 changes: 79 additions & 6 deletions app/src/main/java/com/openathena/DEMParser.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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);
}
}

Expand Down
46 changes: 0 additions & 46 deletions app/src/main/java/com/openathena/DTEDParser.java
View file

This file was deleted.

5 changes: 5 additions & 0 deletions app/src/main/java/com/openathena/MainActivity.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down
29 changes: 14 additions & 15 deletions app/src/main/java/com/openathena/TargetGetter.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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;
Expand Down

0 comments on commit 699befd

Please sign in to comment.