-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathMaxarDtmReader.java
More file actions
1594 lines (1337 loc) · 65.4 KB
/
MaxarDtmReader.java
File metadata and controls
1594 lines (1337 loc) · 65.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// ReadMaxarDtm.java
// Bobby Krupczak
// with ChatGPT help
// javac -cp "lib/*" MaxarDtmReader.java
// java -cp ".:lib/*" MaxarDtmReader
// MaxarDtmReader.java – robust DTM reader using mil.nga.tiff that will
// try a variety of approaches for finding relevant metadata
// because it also will find/use gdalinfo if present
// maxar DTMs encode some params in gdalinfo, OT DTMs in tie points, scales, etc.
// maxar does encode vertical dataum in gdal metadata
// https://chatgpt.com/c/68ea6565-1f40-8330-9669-c5f675a75528
//
// EPSG: 4269 is NAD83
// EPSG: 4326 is WGS84
// NAVD88 [EPSG: 5703] is orthometric height H or AMSL
// WGS84 = EGM96 + offset
//
// h = ellipsoidal height
// H = orthometric height
// N = geoid height or offset
// h = H + N
// wgs84alt = egmalt + offset
//
// EPSG:3855 is EGM2008 or orthometric or AMSL height, COP30, 3DEP (NAVD88), EUDTM,
// EPSG:5773 is EGM96 or orthometric or AMSL height, SRTM, DTED
// EPSG:4979 is WGS84 ellipsoidal height,
// EPSG:5703 is NAVD88 orthometric height
// see this webpage for various WGS84 variants
// https://support.esri.com/en-us/knowledge-base/wgs-1984-is-not-what-you-think-000036058
import mil.nga.tiff.FileDirectory;
import mil.nga.tiff.FileDirectoryEntry;
import mil.nga.tiff.TiffReader;
import mil.nga.tiff.Rasters;
import mil.nga.tiff.TIFFImage;
import mil.nga.tiff.util.TiffException;
import com.agilesrc.dem4j.dted.DTEDLevelEnum;
import com.agilesrc.dem4j.dted.impl.FileBasedDTED;
import com.agilesrc.dem4j.exceptions.CorruptTerrainException;
import com.agilesrc.dem4j.exceptions.InvalidValueException;
import com.agilesrc.dem4j.Point;
import org.locationtech.proj4j.CRSFactory;
import org.locationtech.proj4j.CoordinateReferenceSystem;
import org.locationtech.proj4j.CoordinateTransform;
import org.locationtech.proj4j.CoordinateTransformFactory;
import org.locationtech.proj4j.ProjCoordinate;
import java.io.File;
import java.nio.file.Path;
import java.io.IOException;
import java.lang.reflect.Method;
import java.nio.charset.StandardCharsets;
import java.util.*;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import java.nio.file.Paths;
import java.util.Arrays;
import com.openathena.core.GeoTiffDataType;
import com.openathena.core.EGMOffsetProvider;
import com.openathena.core.EGM96OffsetAdapter;
import com.openathena.core.EGM2008OffsetAdapter;
import com.openathena.core.OpenAthenaCore;
import com.openathena.core.RequestedValueOOBException;
import com.openathena.core.MathUtils;
public class MaxarDtmReader implements AutoCloseable
{
// GeoKey IDs
private static final int KEY_GeographicTypeGeoKey = 2048;
private static final int KEY_ProjectedCSTypeGeoKey = 3072;
private static final int TAG_ModelPixelScale = 33550;
private static final int TAG_ModelTiepoint = 33922;
private static final int TAG_ModelTransformation = 34264;
private static final int TAG_GDAL_METADATA = 42112;
private static final int TAG_GDAL_NODATA = 42113;
private static final double idwPower = 1.875d;
private FileDirectory dir;
private Rasters rasters;
private int width, height;
// Affine: x = a0 + a1*col + a2*row; y = b0 + b1*col + b2*row
private double a0, a1, a2, b0, b1, b2;
private double inv00, inv01, inv10, inv11;
private boolean georeferenced;
private Double noData;
private String dataEpsg;
private CoordinateReferenceSystem dataCRS, wgs84;
private CoordinateTransform dataToWgs, wgsToData;
private String gdal; // gdalinfo field if present
private Method mGetPixelSampleDouble; // (x,y,band)->double
private Method mGetFirstPixelSample; // (x,y)->double
private Method mGetPixelSample; // (x,y,band)->Number
private Method mGetPixel; // (x,y)->Object (e.g., Number[] or array)
private String horizontalCRS; // e.g., "EPSG:32616"
private String verticalCRS; // e.g., "ELLIPSOID (WGS84 G1674)" or "EPSG:5703"
private String verticalDatum; // set in testVerticalDatum()
// PixelIsPoint? true = centers, false = corners
private boolean centerAnchored = false;
private final CRSFactory crsFactory = new CRSFactory();
private final CoordinateTransformFactory ctf = new CoordinateTransformFactory();
// Tiny cache so we don’t rebuild transforms for repeated calls
private Map<String, CoordinateTransform> reprojCache = new java.util.HashMap<>();
// once we've picked a directory, keep map of its tag to entry mapping
private Map<Integer,FileDirectoryEntry> directoryIndexTags = new HashMap<Integer,FileDirectoryEntry>();
// corners of the bounding box calculated from the GeoTiff or DTED, not
// from the filename; lat,lon
private double n,s,e,w;
protected EGMOffsetProvider offsetProvider = null; // initialized later
// DTED params
public int numRows, numCols;
public double latSpacing, lonSpacing;
// public params
protected transient File geofile;
private transient TIFFImage tiff;
private FileBasedDTED dted;
public boolean isDTED = false;
public String filepath;
public String filename;
public String extension;
public DTEDLevelEnum dtedLevel = null;
public int numAltLookups = 0;
// type of GeoTiff this object represents
public GeoTiffDataType gType;
public MaxarDtmReader(File geofile) throws Exception
{
this.geofile = geofile;
filepath = geofile.getPath();
filename = geofile.toPath().getFileName().toString();
int dotIndex = filepath.lastIndexOf('.');
extension = filepath.substring(dotIndex +1);
gType = GeoTiffDataType.fromExtension(extension);
this.geofile = geofile;
// check if it is a DTED
if (looksLikeDTED(geofile.toPath())) {
// System.out.println(filename + " is DTED");
isDTED = true;
// read DTED instead and return;
readDted(filepath);
}
else {
// if maxar, we'll discover it by finding gdalinfo in tiff
// and set the gType even if extension is .tif[f]
readGeofile(filepath);
}
}
// read a GeoTiff file
private void readGeofile(String filepath) throws Exception
{
// System.out.println("readGeofile: "+filepath);
long t0, t1;
tiff = TiffReader.readTiff(geofile);
// Gather all IFDs (root + subIFDs if exposed)
List<FileDirectory> all = collectAllDirectories(tiff);
// System.out.println("readGeoFile: number of directories is "+all.size());
// could be multiple directories; examine and pick best (that has georeference)
FileDirectory picked = null;
long bestArea = -1;
boolean bestHas = false;
for (FileDirectory d : all) {
int w = d.getImageWidth().intValue();
int h = d.getImageHeight().intValue();
//System.out.println("Examining directory "+d);
long area = (long) w * (long) h;
boolean has = hasAnyGeorefSignature(d);
if (picked == null || (has && !bestHas) || (has == bestHas && area > bestArea)) {
picked = d; bestArea = area; bestHas = has;
// System.out.println("Picked best");
break; // no need to continue
}
}
if (picked == null) picked = tiff.getFileDirectory();
this.dir = picked;
this.centerAnchored = isCenterAnchored(dir);
this.width = dir.getImageWidth().intValue();
this.height = dir.getImageHeight().intValue();
// Build affine params for mapping between x,y and pixel array
AffineParams ap = tryBuildAffineRobust(dir);
if (ap != null) {
setAffine(ap.a0, ap.a1, ap.a2, ap.b0, ap.b1, ap.b2);
this.georeferenced = true;
} else {
setAffine(0,1,0, 0,0,-1);
this.georeferenced = false;
System.err.println("[MaxarDtmReader] No georeferencing found; operating in pixel space.");
}
// this takes longer for cop30 and not sure why
this.rasters = dir.readRasters();
// Horizontal & vertical CRS/datum detection
this.horizontalCRS = determineHorizontalCRS(dir);
this.verticalCRS = determineVerticalCRS(dir);
if (georeferenced) {
// Prefer the horizontal CRS we parsed earlier
this.dataEpsg = (this.horizontalCRS != null) ? this.horizontalCRS : detectEpsgRobust(dir);
if (this.dataEpsg != null) {
// this will create transforms
enableCrs(this.dataEpsg);
} else {
System.err.println("[MaxarDtmReader] Warning: horiz CRS not found; will initialize transforms lazily on first query.");
}
}
// determine corners n,s,e,w
Bounds wgs = getBoundsWGS84();
if (wgs != null) {
this.w = wgs.minX;
this.e = wgs.maxX;
this.s = wgs.minY;
this.n = wgs.maxY;
}
// see what file said was vertical datum and potentially override
// the gType vertical datum; if needed, set offset provider
// takes a long time on 3dep, cop30, eudtm
// is this due to size of EGM2008? Yes
testVerticalDatum();
} // readGeofile
// read/process DTED file and get its parameters
// only handles DTEDs
private void readDted(String inputFilePath)
{
// System.out.println("readDted: "+filepath);
try {
File geofile = new File(inputFilePath);
dted = new FileBasedDTED(geofile);
dtedLevel = dted.getDTEDLevel();
// System.out.println("DTED level "+dtedLevel);
if (dtedLevel.equals(DTEDLevelEnum.DTED0) || dtedLevel.equals(DTEDLevelEnum.DTED1)) {
System.out.println("DTED2 or DTED3 or higher is required");
throw new TiffException("DTED2, 3 or higher is required");
}
// can always test this against gdalinfo output
n = dted.getNorthWestCorner().getLatitude();
w = dted.getNorthWestCorner().getLongitude();
s = dted.getSouthEastCorner().getLatitude();
e = dted.getSouthEastCorner().getLongitude();
latSpacing = dted.getLatitudeInterval();
lonSpacing = dted.getLongitudeInterval();
numRows = dted.getRows();
numCols = dted.getColumns();
//System.out.println("dted: n: "+n);
//System.out.println("dted: w: "+w);
//System.out.println("dted: s: "+s);
//System.out.println("dted: e: "+e);
//System.out.println("dted: resolution "+dted.getResolution().getRows()+","
// +dted.getResolution().getColumns()+" "+dted.getResolution().getSpacing()+" degrees");
this.isDTED = true;
this.gType = GeoTiffDataType.DTED2;
testVerticalDatum();
}
catch (Exception e) {
throw new TiffException(e.getMessage());
}
return;
} // readDted
// once DEM has been parsed and loaded, check if vertical CRS was dectected in metadata
// and possiblyl override default that was based on file extension
// Maxar DTMs have vertical CRS set in GDAL info
private void testVerticalDatum()
{
// default to gType
offsetProvider = gType.getOffsetProvider();
verticalDatum = gType.getVertDatum();
// look at verticalCRS string and if set, override gType
if ((verticalCRS == null) || (verticalCRS.equals(""))) {
return;
}
if (verticalCRS.equals("EPSG:0")) {
return;
}
// System.out.println("testVerticalDatum: overriding vertical datum with "+verticalCRS);
switch (verticalCRS) {
case "EPSG:3855":
offsetProvider = new EGM2008OffsetAdapter();
verticalDatum = verticalCRS;
break;
case "EPSG:5773":
offsetProvider = new EGM96OffsetAdapter();
verticalDatum = verticalCRS;
break;
case "EPSG:5703":
// NAVD88 -> EGM2008 is acceptable substitute
offsetProvider = new EGM2008OffsetAdapter();
verticalDatum = verticalCRS;
break;
case "EPSG:4979":
// WGS84; no need for offset provider
// but we load 2008 so we can also provide egm2008 altitude
offsetProvider = new EGM2008OffsetAdapter();
verticalDatum = verticalCRS;
break;
default:
System.out.println("Unrecognized vertical datum "+verticalCRS);
}
} // testVerticalDatum
// public api functions
public String getHorizontalCRS() { return horizontalCRS; }
public String getVerticalCRS() { return verticalCRS; }
public int getWidth() { return width; }
public int getHeight() { return height; }
public boolean isGeoreferenced() { return georeferenced; }
public String getDataEpsg() { return dataEpsg; }
// public Optional<Double> getNoData() { return Optional.ofNullable(noData); }
public String getVerticalDatum() { return verticalDatum; }
public String getHorizontalDatum() { return gType.getHorizDatum(); }
public String getFilename() { return filename; }
public String getFilepath() { return filepath; }
public String getExtension() { return extension; }
public double getS() { return s; }
public double getMinLat() { return s; }
public double getE() { return e; }
public double getMaxLon() { return e; }
public double getW() { return w; }
public double getMinLon() { return w; }
public double getN() { return n; }
public double getMaxLat() { return n; }
// Bounds in the raster's native CRS using a center-aware corner convention
public Bounds getBoundsDataCRS()
{
requireGeoref("Bounds in data CRS");
// Decide pixel coordinates for the four corners based on anchoring
final boolean center = this.centerAnchored; // set this in ctor from isCenterAnchored(dir)
final double x0 = center ? -0.5 : 0.0;
final double y0 = center ? -0.5 : 0.0;
final double x1 = center ? (width - 0.5) : width;
final double y1 = center ? (height - 0.5) : height;
double[] p00 = worldFromPixel(x0, y0); // UL
double[] p10 = worldFromPixel(x1, y0); // UR
double[] p01 = worldFromPixel(x0, y1); // LL
double[] p11 = worldFromPixel(x1, y1); // LR
double minX = Math.min(Math.min(p00[0], p10[0]), Math.min(p01[0], p11[0]));
double maxX = Math.max(Math.max(p00[0], p10[0]), Math.max(p01[0], p11[0]));
double minY = Math.min(Math.min(p00[1], p10[1]), Math.min(p01[1], p11[1]));
double maxY = Math.max(Math.max(p00[1], p10[1]), Math.max(p01[1], p11[1]));
return new Bounds(minX, minY, maxX, maxY);
}
public Bounds getBoundsWGS84Old() {
requireGeoref("Bounds in WGS84");
Bounds b = getBoundsDataCRS();
ProjCoordinate p1 = transform(dataToWgs, b.minX, b.minY);
ProjCoordinate p2 = transform(dataToWgs, b.minX, b.maxY);
ProjCoordinate p3 = transform(dataToWgs, b.maxX, b.minY);
ProjCoordinate p4 = transform(dataToWgs, b.maxX, b.maxY);
double minLon = Math.min(Math.min(p1.x, p2.x), Math.min(p3.x, p4.x));
double maxLon = Math.max(Math.max(p1.x, p2.x), Math.max(p3.x, p4.x));
double minLat = Math.min(Math.min(p1.y, p2.y), Math.min(p3.y, p4.y));
double maxLat = Math.max(Math.max(p1.y, p2.y), Math.max(p3.y, p4.y));
return new Bounds(minLon, minLat, maxLon, maxLat);
}
// get bounds in WGS84 (lon/lat) decimal degrees of this DEM
public Bounds getBoundsWGS84()
{
requireGeoref("Bounds in WGS84");
Bounds b = getBoundsDataCRS();
org.locationtech.proj4j.ProjCoordinate src = new org.locationtech.proj4j.ProjCoordinate();
org.locationtech.proj4j.ProjCoordinate dst = new org.locationtech.proj4j.ProjCoordinate();
// Transform the same four corners used above
double[][] xy = {
{ b.minX, b.minY }, { b.minX, b.maxY }, { b.maxX, b.minY }, { b.maxX, b.maxY }
};
double west = Double.POSITIVE_INFINITY, east = Double.NEGATIVE_INFINITY;
double south = Double.POSITIVE_INFINITY, north = Double.NEGATIVE_INFINITY;
for (double[] p : xy) {
src.x = p[0]; src.y = p[1];
dataToWgs.transform(src, dst); // lon,lat
west = Math.min(west, dst.x);
east = Math.max(east, dst.x);
south = Math.min(south, dst.y);
north = Math.max(north, dst.y);
}
return new Bounds(west, south, east, north);
}
// use bilinear interpolation to get elevation at x,y
private double getElevationProjectedBilinear(double x, double y)
{
double[] rc = pixelFromWorld(x, y);
double col = rc[0], row = rc[1];
if (col < 0 || row < 0 || col > (width-1) || row > (height-1)) return Double.NaN;
int c0 = (int)Math.floor(col), r0 = (int)Math.floor(row);
int c1 = Math.min(c0+1, width-1), r1 = Math.min(r0+1, height-1);
double dc = col - c0, dr = row - r0;
Double z00 = sample(c0, r0), z10 = sample(c1, r0), z01 = sample(c0, r1), z11 = sample(c1, r1);
if (z00==null || z10==null || z01==null || z11==null) return Double.NaN;
double z0 = z00*(1-dc) + z10*dc;
double z1 = z01*(1-dc) + z11*dc;
return z0*(1-dr) + z1*dr;
}
// Convenience: IDW with radius=1 (3x3 window), power=idwPower, epsilon=1e-12.
// for use with geotiffs only
private double getElevationProjectedIDW(double xData, double yData)
{
return getElevationProjectedIDW(xData, yData, /*radius*/ 1, idwPower, /*epsilon*/ 1e-12);
}
/**
* Inverse distance weighting (IDW) interpolation in the raster’s native CRS.
* for use with geotiffs only
*
* @param xData X/easting in the DTM's data CRS
* @param yData Y/northing in the DTM's data CRS
* @param radius neighborhood radius in pixels (1 = 3x3, 2 = 5x5, etc.)
* @param power IDW power parameter (common: 1.5–3; 2 is typical)
* @param epsilon small distance threshold; if the query is essentially on a pixel center,
* return that pixel’s value directly (avoids division by ~0).
*/
private double getElevationProjectedIDW(double xData, double yData, int radius, double power, double epsilon)
{
requireGeoref("IDW elevation");
// Convert world → pixel coordinates (col,row), fractional.
double[] rc = pixelFromWorld(xData, yData);
double col = rc[0], row = rc[1];
// If OOB (outside pixel-edge box [0..width]x[0..height]), short-circuit
if (col < 0 || row < 0 || col > width || row > height) {
return Double.NaN;
}
// Center-aware distances: pixel centers are at (c+0.5, r+0.5) when your affine is corner-based.
int c0 = (int) Math.floor(col);
int r0 = (int) Math.floor(row);
int cMin = Math.max(0, c0 - radius);
int cMax = Math.min(width - 1, c0 + radius);
int rMin = Math.max(0, r0 - radius);
int rMax = Math.min(height - 1, r0 + radius);
double wsum = 0.0;
double vsum = 0.0;
for (int r = rMin; r <= rMax; r++) {
for (int c = cMin; c <= cMax; c++) {
// Fetch the pixel; assume your existing 'sample(c,r)' returns null for NoData/OOB
Double v = sample(c, r);
if (v == null) continue;
// Distance in pixel space to the pixel center:
double dc = (col - (c + 0.5));
double dr = (row - (r + 0.5));
double d2 = dc*dc + dr*dr;
// Exact hit (or virtually exact): return the pixel value immediately
if (d2 <= epsilon * epsilon) {
return v;
}
double d = Math.sqrt(d2);
double w = 1.0 / Math.pow(d, power);
wsum += w;
vsum += w * v;
}
}
if (wsum == 0.0) return Double.NaN;
return vsum / wsum;
} // getElevationProjectedIDW
// for a lat,lon, use our EGM96 offset provider to return the lat,lon
// we should allocate an offsetprovider anyway ?
// Maxar DEMs don't need offset provider but we load EGM2008 for conversion
// purposes
public double getEGMOffsetForLatLon(double lat, double lon)
{
return offsetProvider.getEGMOffsetAtLatLon(lat,lon);
}
// we could be smarter here because many of our altitudes are
// in EGM but its too easy to just get the altitude in WGS84 and subtract
// the offset; less code to maintain but potentially calls getEGMOffset twice
// maxar DEMs don't need an offset provider since their vertical datum is wgs84
// but we load EGM2008 for maxar so we can convert altitudes
public double getEGMAltFromLatLon(double lat, double lon) throws RequestedValueOOBException, CorruptTerrainException
{
double wgs84alt = getAltFromLatLon(lat,lon);
double offset = getEGMOffsetForLatLon(lat,lon);
return wgs84alt - offset;
}
// given a DTED, and lat,lon in decimal degrees, return WGS84 altitude
// no interpolation; agilesrc dem4j takes the nearest point and returns that elevation
private double getAltFromLatLonDted(double lat, double lon) throws RequestedValueOOBException, CorruptTerrainException
{
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
// re issue #54
double WGS84_altitude = EGM96_altitude + offsetProvider.getEGMOffsetAtLatLon(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);
}
}
// use inverse distance weight with X neighbors/elevations to calculate altitude
// if it turns out that the targetLat,targetLon is exact, return that value instead
// of IDW
private double idwInterpolation(double targetLat, double targetLon, Point[] neighbors, double[] elevations, double power)
{
double sumWeights = 0.0;
double sumWeightedElevations = 0.0;
int i;
// System.out.println("idwInterpolation: target is "+targetLat+","+targetLon);
// System.out.println("idwInterpolation: neighbors are "+neighbors);
// System.out.println("idwInterpolation: elevations are "+Arrays.toString(elevations));
for (i=0; i<neighbors.length; i++) {
// System.out.println("idwInterpolation: neighbor "+i+" "+neighbors[i].getLatitude()+","+neighbors[i].getLongitude());
double distance = MathUtils.haversine(targetLon, targetLat, neighbors[i].getLongitude(), neighbors[i].getLatitude(), elevations[i]);
double weight = 1.0d / Math.pow(distance, power);
// System.out.println("idwInterpolation: distance: "+distance+" weight: "+weight);
// if distance is ~= 0.0 then we got lucky and the point is actually target;
// if so, return that alt w/o need to interpolate
if (distance <= 0.1) {
return elevations[i];
}
sumWeights += weight;
sumWeightedElevations += weight * elevations[i];
}
// System.out.println("idwInterpolation: returning "+sumWeightedElevations / sumWeights);
return sumWeightedElevations / sumWeights;
}
// get WGS84 altitude from lat,lon, via DTED, using IDW
private double getAltFromLatLonDtedIDW(double lat, double lon) throws RequestedValueOOBException, CorruptTerrainException
{
// target point
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()) {
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
// re issue #180, fix incorrect equation for applying geoid offset
double WGS84_altitude = EGM96_altitude + offsetProvider.getEGMOffsetAtLatLon(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 {
// 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();
// Perform IDW interpolation over points 1-4 and elevations 1-4
// code taken from Android OpenAthena and ported here
Point[] neighbors = new Point[] { p1, p2, p3, p4};
double[] elevations = new double[] { e1, e2, e3, e4};
double interpolatedAltitude = idwInterpolation(lat,lon,neighbors,elevations,idwPower);
// System.out.println("interpolatedAltitude is "+interpolatedAltitude);
// System.out.println("non interpolated Altitude is "+ dted.getElevation(point).getElevation());
// Convert from EGM96 AMSL orthometric height to WGS84 HAE
// re issue #180, fix incorrect equation for applying geoid offset
double WGS84_altitude = interpolatedAltitude + offsetProvider.getEGMOffsetAtLatLon(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);
}
}
// main public API for getting elevation for lat,lon
// get WGS84 altitude for lat,lon decimal degrees; calls for
// both DTED and Geotiff;
// uses appropriate IDW with surrounding neighbors to calculate
// altitude
public double getAltFromLatLon(double latDeg, double lonDeg) throws RequestedValueOOBException, CorruptTerrainException
{
if (this.isDTED) {
return getAltFromLatLonDtedIDW(latDeg, lonDeg);
}
requireGeoref("Elevation (lat,lon in WGS84)");
// Ensure we know/enable the data CRS
if (this.dataCRS == null) {
// Try, in order: previously detected EPSG, parsed horizontalCRS, or parse now
String epsg = (this.dataEpsg != null) ? this.dataEpsg
: (this.horizontalCRS != null ? this.horizontalCRS : determineHorizontalCRS(this.dir));
if (epsg == null || epsg.isBlank()) {
throw new IllegalStateException("Data CRS unknown; cannot reproject from WGS84.");
}
enableCrs(epsg); // sets dataCRS, wgs84, dataToWgs, wgsToData, dataEpsg
}
// proj4j expects (lon, lat) when transforming geographic coords
ProjCoordinate p = transform(this.wgsToData, lonDeg, latDeg);
// double alt = getElevationProjectedBilinear(p.x, p.y);
double alt = getElevationProjectedIDW(p.x,p.y);
if (Double.isNaN(alt)) {
throw new RequestedValueOOBException("getAltFromLatLon args out of bounds due to min/max lat/lon!",latDeg,lonDeg);
}
// return must be in WGS84 HAE
switch (verticalDatum) {
case "EPSG:3855": // EGM2008
case "EPSG:5773": // EGM96
case "EPSG:5703": // NAVD88 -> approx with EGM2008
// offset provider better be set!!
double offset = offsetProvider.getEGMOffsetAtLatLon(latDeg,lonDeg);
alt = alt + offset;
case "EPSG:4979":
// do nothing as vertical datum is WGS84 already
}
return alt;
}
@Override public void close() {}
// internal methods for parsing/manipulating GeoTiff
private static double[] readDoubleArrayTag(Map<Integer,FileDirectoryEntry> t, int tagId)
{
FileDirectoryEntry e = t.get(tagId);
if (e == null) return null;
Object v = e.getValues();
if (v == null) return null;
//System.out.println("readDoubleArrayTag: read object");
//System.out.println("readDoubleArrayTag: object is of type "+v.getClass().getName());
List<?> L = (List<?>) v;
double[] out = new double[L.size()];
for (int i = 0; i < L.size(); i++) {
Object o = L.get(i);
if (o == null) { out[i] = Double.NaN; continue; }
if (o instanceof Number) {
out[i] = ((Number) o).doubleValue();
} else {
// fallback: parse numeric strings
String s = o.toString().trim();
try {
out[i] = Double.parseDouble(s);
} catch (NumberFormatException nfe) {
out[i] = Double.NaN; // or throw, if you prefer strict behavior
}
}
}
return out;
}
private void setAffine(double A0,double A1,double A2,double B0,double B1,double B2)
{
this.a0=A0; this.a1=A1; this.a2=A2; this.b0=B0; this.b1=B1; this.b2=B2;
double det = a1*b2 - a2*b1;
if (Math.abs(det) < 1e-12) throw new IllegalStateException("Non-invertible geotransform (det≈0)");
this.inv00 = b2/det; this.inv01 = -a2/det;
this.inv10 = -b1/det; this.inv11 = a1/det;
}
private void enableCrs(String epsg)
{
CRSFactory cf = new CRSFactory();
this.dataCRS = cf.createFromName(epsg);
this.wgs84 = cf.createFromName("EPSG:4326");
CoordinateTransformFactory ctf = new CoordinateTransformFactory();
this.dataToWgs = ctf.createTransform(dataCRS, wgs84);
this.wgsToData = ctf.createTransform(wgs84, dataCRS);
this.dataEpsg = epsg;
}
private void requireGeoref(String op) {
if (!georeferenced) throw new IllegalStateException(op + " requires georeferencing.");
}
private static List<FileDirectory> collectAllDirectories(TIFFImage tiffImage)
{
return tiffImage.getFileDirectories();
}
// PixelIsPoint ⇒ center-anchored; PixelIsArea ⇒ corner-anchored. Default: Area (false).
private boolean isCenterAnchored(FileDirectory d)
{
if (directoryIndexTags.size() == 0) {
indexDirectoryTags(d);
}
// Prefer GDAL metadata if present
gdal = findGdalMetadata(d);
if (gdal != null) {
gType = GeoTiffDataType.MAXAR;
String v = findMetadataItem(gdal, "AREA_OR_POINT");
if (v != null) return v.trim().equalsIgnoreCase("Point");
}
// Fall back to scanning ASCII entries
for (FileDirectoryEntry e : safeEntries(d)) {
String s = toAscii(e.getValues());
if (s == null) continue;
String u = s.toUpperCase();
if (u.contains("AREA_OR_POINT=POINT")) return true;
if (u.contains("AREA_OR_POINT=AREA")) return false;
}
// Heuristic: a single tiepoint at (i,j) ≈ (0.5,0.5) is also a strong signal for center
// double[] tp = tryGetModelTiepoint(d);
double[] tp = readDoubleArrayTag(directoryIndexTags, TAG_ModelTiepoint);
if (tp != null && tp.length >= 2) {
double i = tp[0], j = tp[1];
if (Math.abs(i - 0.5) < 1e-9 && Math.abs(j - 0.5) < 1e-9) return true;
}
return false;
}
// does this tiff/dir have any georef signatures
// srtm,cop30,3dep have model pixel scale and tiepoints
// max has double16matrix and gdal metadata
private boolean hasAnyGeorefSignature(FileDirectory d)
{
if (directoryIndexTags.size() == 0) {
indexDirectoryTags(d);
}
if ((readDoubleArrayTag(directoryIndexTags, TAG_ModelPixelScale) != null) &&
(readDoubleArrayTag(directoryIndexTags, TAG_ModelTiepoint) != null)) {
// System.out.println("hasAnyGeorefSignature: model tie point and pixel scale");
return true;
}
if (findDouble16Matrix(d) != null) {
// System.out.println("hasAnyGeorefSignature: double16 matrix");
return true;
}
if (extractGeoTransformFromGdalMetadata(d) != null) {
// System.out.println("hasAnyGeorefSignature: gdal metadata");
return true;
}
//System.out.println("hasAnyGeorefSignature: no");
return false;
}
// small helper near your other utils
private static boolean near(double v, double tgt) { return Math.abs(v - tgt) < 1e-9; }
/** Returns corner-based affine (GDAL convention) from Scale+Tiepoint.
* Applies half-pixel shift when tiepoint is at (0.5,0.5) or when AREA_OR_POINT=Point.
*/
private AffineParams affineFromScaleTiepointCornered(double[] scale, double[] tie, FileDirectory d)
{
// scale = [sx, sy, (sz?)]; tie = [i, j, k, x, y, z] (use the first tiepoint)
double sx = scale[0], sy = scale[1];
double i = tie[0], j = tie[1], x = tie[3], y = tie[4];
// Start from generic formulas
double a1 = sx, a2 = 0.0;
double b1 = 0.0, b2 = -sy;
// Is this center-anchored?
boolean centerTP = Math.abs(i - 0.5) < 1e-9 && Math.abs(j - 0.5) < 1e-9;
boolean pixelIsPoint = isPixelIsPoint(d); // checks GDAL_METADATA or TIFF AREA_OR_POINT
boolean needsCenterToCorner = centerTP || pixelIsPoint;
double a0, b0;
if (needsCenterToCorner) {
// Convert center origin → corner origin
a0 = (x - 0.5 * sx) - (i - 0.5) * sx;
b0 = (y + 0.5 * sy) + (j - 0.5) * sy;
} else {
// Assume tiepoint already references a pixel corner
a0 = x - i * sx;
b0 = y + j * sy;
}
return new AffineParams(a0, a1, a2, b0, b1, b2);
}
// Detect PixelIsPoint (true) vs PixelIsArea (false). Defaults to false if unknown.
private boolean isPixelIsPoint(FileDirectory d) {
// 1) Look in GDAL_METADATA <Item name="AREA_OR_POINT">Point|Area</Item>
gdal = findGdalMetadata(d);
if (gdal != null) {
gType = GeoTiffDataType.MAXAR;
String v = findMetadataItem(gdal, "AREA_OR_POINT");
if (v != null) return v.trim().equalsIgnoreCase("Point");
}
// 2) Some datasets set TIFF tag 'AREA_OR_POINT' as ASCII elsewhere; scan entries
for (FileDirectoryEntry e : safeEntries(d)) {
String s = toAscii(e.getValues());
if (s != null && s.toUpperCase().contains("AREA_OR_POINT=POINT")) return true;
if (s != null && s.toUpperCase().contains("AREA_OR_POINT=AREA")) return false;
}
return false; // default: treat as PixelIsArea (corner)
}
/** Corner-based affine from a 4x4 ModelTransformation matrix.
* Many writers map pixel centers; apply half-pixel shift to match GDAL corners.
*/
private static AffineParams affineFromMatrixCornered(double[] m16)
{
// Row-major 4x4: we only need the 2D part
double a1 = m16[0], a2 = m16[1], a0 = m16[3];
double b1 = m16[4], b2 = m16[5], b0 = m16[7];
// Shift center → corner
double adjA0 = a0 - 0.5 * a1 - 0.5 * a2;
double adjB0 = b0 - 0.5 * b1 - 0.5 * b2;
return new AffineParams(adjA0, a1, a2, adjB0, b1, b2);
}
/**
* Build a corner-based (GDAL-style) affine:
* Order of preference:
* 1) GDAL_METADATA <GeoTransform> (already corner-based)
* 2) ModelPixelScale + ModelTiepoint (center→corner if (0.5,0.5) or PixelIsPoint)
* 3) ModelTransformation (4×4) (center→corner shift)
* 4) Scans for 16-dbl matrix or scale/tie arrays with same normalization
*/
private AffineParams tryBuildAffineRobust(FileDirectory d)
{
if (directoryIndexTags.size() == 0) {
indexDirectoryTags(d);
}
// 1) Preferred: GDAL <GeoTransform> (exactly matches gdalinfo)
double[] gt = extractGeoTransformFromGdalMetadata(d); // your existing helper
if (gt != null && gt.length == 6) {
// System.out.println("tryBuildAffineRobust: using geoTransform from gdal");
return new AffineParams(gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
}
// 2) Scale + Tiepoint (normalize to corner)
//double[] scale = tryGetModelPixelScale(d); // your existing helper
//double[] tie = tryGetModelTiepoint(d); // your existing helper
double[] scale = readDoubleArrayTag(directoryIndexTags, TAG_ModelPixelScale);
double[] tie = readDoubleArrayTag(directoryIndexTags, TAG_ModelTiepoint);
if (scale != null && scale.length >= 2 && tie != null && tie.length >= 6) {
double sx = scale[0], sy = scale[1];
double i = tie[0], j = tie[1], x = tie[3], y = tie[4];
// Start from indices → world
double originX = x - i * sx;
double originY = y + j * sy;
// If center-anchored (either (0.5,0.5) OR PixelIsPoint), shift to corner
boolean centerIJ = near(i, 0.5) && near(j, 0.5);
if (centerIJ || isPixelIsPoint(d)) {
originX -= 0.5 * sx;
originY += 0.5 * sy; // note: we'll set b2 negative below
}
// North-up (no rotation) assumption for Scale+Tiepoint
// System.out.println("tryBuildAffineRobust: using scale and tiepoints");
return new AffineParams(originX, sx, 0.0, originY, 0.0, -sy);
}
// 3) ModelTransformation (4×4) → shift center→corner
// we don't have getModelTransform method in our mil.nga.tiff
// 4a) Scan for another 16-dbl matrix and normalize the same way
double[] mt2 = findDouble16Matrix(d); // your existing helper
if (mt2 != null && mt2.length == 16) {
//System.out.println("tryBuildAffineRobust: using double 16 matrix");