如何使用 java 从 GeoTIFF 图像中提取纬度和经度值?

Question

我有一个卫星图像，我必须处理它。到目前为止，我已经使用ImageIO和JAI读取了该文件。我现在希望找到与存储在TIFF文件中的每个像素相对应的纬度和经度值。提前谢谢。

Answer 1

假设您有一个GeoTiff，即具有地理参考的GeoTiff，那么GeoTools可以为您处理它。

下面显示了如何使用实际坐标访问像素：

AbstractGridFormat format = GridFormatFinder.findFormat(parameters);
AbstractGridCoverage2DReader reader = format.getReader(parameters);
System.out.println(parameters);
CoordinateReferenceSystem crs = reader.getCoordinateReferenceSystem();
System.out.println(crs);
GridCoverage2D cov = null;
try {
cov = reader.read(null);
System.out.println("type: " + cov.getName());
Envelope2D bbox = cov.getEnvelope2D();
me.printBBOX(bbox);
int nOver = cov.getNumOverviews();
if (nOver > 0) {
System.out.println("" + nOver + " overviews");
}
GridEnvelope gridRange2D = cov.getGridGeometry().getGridRange();
System.out.println("dimension: " + gridRange2D.getDimension());
for (int i = 0; i < gridRange2D.getDimension(); i++) {
System.out
.println("dimension " + i + " pixel range " + gridRange2D.getLow(i) + " - " + gridRange2D.getHigh(i));
}
CoordinateReferenceSystem wgs84 = DefaultGeographicCRS.WGS84;
CoordinateReferenceSystem target = cov.getCoordinateReferenceSystem();// CRS.decode("EPSG:3997",
    // true);
GeometryFactory gf = new GeometryFactory();
MathTransform targetToWgs = CRS.findMathTransform(target, wgs84);
int count = 0;
for (int j = gridRange2D.getLow(1); j < gridRange2D.getHigh(1); j++) {
for (int i = gridRange2D.getLow(0); i < gridRange2D.getHigh(0); i++) {
if (count++ > 10)
System.exit(1);
GridCoordinates2D coord = new GridCoordinates2D(i, j);
DirectPosition p = cov.getGridGeometry().gridToWorld(coord);
Point point = gf.createPoint(new Coordinate(p.getOrdinate(0), p.getOrdinate(1)));
Geometry wgsP = JTS.transform(point, targetToWgs);
System.out.format("(%d %d) -> POINT(%.2f %.2f) -> POINT(%.2f %.2f)%n", i, j, point.getCoordinate().x,
point.getCoordinate().y, wgsP.getCentroid().getCoordinate().x, wgsP.getCentroid().getCoordinate().y);
}
}
} catch (IOException giveUp) {
throw new RuntimeException(giveUp);
} catch (FactoryException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (TransformException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
} catch (ParseException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}

如果你只需要边界框，那么你可以使用这样的东西：

private void printBBOX(Envelope2D bbox) {
MathTransform trans;
ReferencedEnvelope boxWGS = null;
CoordinateReferenceSystem crs = bbox.getCoordinateReferenceSystem();
try {
trans = CRS.findMathTransform(crs, DefaultGeographicCRS.WGS84);
boxWGS = (ReferencedEnvelope) JTS.toGeographic(JTS.toEnvelope(bbox), crs);
} catch (FactoryException e) {
// TODO Auto-generated catch block
e.printStackTrace();
return;
} catch (TransformException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
CoordinateFormat formatter = new CoordinateFormat();
formatter.setNumberPattern("####0.00");
formatter.setAnglePattern("DD.dd");
formatter.setCoordinateReferenceSystem(DefaultGeographicCRS.WGS84);
Format xFormat = formatter.getFormat(0);
Format yFormat = formatter.getFormat(1);
String unit = "" + DefaultGeographicCRS.WGS84.getAxis(0).getUnit().toString();
int x = 0;
int y = 1;
int wx = 0;
int wy = 1;
char[][] label = { { 'S', 'N' }, { 'W', 'E' } };
if (CRS.getAxisOrder(crs) == AxisOrder.NORTH_EAST) {
x = 1;
y = 0;
}
if (CRS.getAxisOrder(DefaultGeographicCRS.WGS84) == AxisOrder.NORTH_EAST) {
wx = 1;
wy = 0;
}
String format = "%s:t( %.4f, %.4f) (%.4f, %.4f) (%s, %s)%n";
System.out.format(format, "Lower Left", bbox.getMinimum(x), bbox.getMinimum(y), boxWGS.getMinimum(wx),
boxWGS.getMinimum(wy),
(xFormat.format(Math.abs(boxWGS.getMinimum(wx))) + unit
+ (boxWGS.getMinimum(wx) < 0 ? label[wx][0] : label[wx][1])),
(yFormat.format(Math.abs(boxWGS.getMinimum(wy))) + unit)
+ (boxWGS.getMinimum(wy) < 0 ? label[wy][0] : label[wy][1]));
System.out.format(format, "Upper Left", bbox.getMinimum(x), bbox.getMaximum(y), boxWGS.getMinimum(wx),
boxWGS.getMaximum(wy),
(xFormat.format(Math.abs(boxWGS.getMinimum(wx))) + unit
+ (boxWGS.getMinimum(wx) < 0 ? label[wx][0] : label[wx][1])),
(yFormat.format(Math.abs(boxWGS.getMaximum(wy))) + unit)
+ (boxWGS.getMaximum(wy) < 0 ? label[wy][0] : label[wy][1]));
System.out.format(format, "Lower Right", bbox.getMaximum(x), bbox.getMinimum(y), boxWGS.getMaximum(wx),
boxWGS.getMinimum(wy),
(xFormat.format(Math.abs(boxWGS.getMaximum(wx))) + unit
+ (boxWGS.getMaximum(wx) < 0 ? label[wx][0] : label[wx][1])),
(yFormat.format(Math.abs(boxWGS.getMinimum(wy))) + unit)
+ (boxWGS.getMinimum(wy) < 0 ? label[wy][0] : label[wy][1]));
System.out.format(format, "Upper Right", bbox.getMaximum(x), bbox.getMaximum(y), boxWGS.getMaximum(wx),
boxWGS.getMaximum(wy),
(xFormat.format(Math.abs(boxWGS.getMaximum(wx))) + unit
+ (boxWGS.getMaximum(wx) < 0 ? label[wx][0] : label[wx][1])),
(yFormat.format(Math.abs(boxWGS.getMaximum(wy))) + unit)
+ (boxWGS.getMaximum(wy) < 0 ? label[wy][0] : label[wy][1]));
}

在我的测试文件中，这会产生以下输出：

//data/natural_earth/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif
GEOGCS["WGS 84", 
DATUM["World Geodetic System 1984", 
SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
AUTHORITY["EPSG","6326"]], 
PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
UNIT["degree", 0.017453292519943295], 
AXIS["Geodetic latitude", NORTH], 
AXIS["Geodetic longitude", EAST], 
AUTHORITY["EPSG","4326"]]
type: geotiff_coverage
Lower Left: ( -90.0000, -180.0000) (-90.0000, -180.0000) (90°S, 180°W)
Upper Left: ( -90.0000, 180.0000) (-90.0000, 180.0000) (90°S, 180°E)
Lower Right:    ( 90.0000, -180.0000) (90.0000, -180.0000) (90°N, 180°W)
Upper Right:    ( 90.0000, 180.0000) (90.0000, 180.0000) (90°N, 180°E)
dimension: 2
dimension 0 pixel range 0 - 21599
dimension 1 pixel range 0 - 10799
(0 0) -> POINT(-179.99 89.99) -> POINT(89.99 -179.99)
(1 0) -> POINT(-179.98 89.99) -> POINT(89.99 -179.98)
(2 0) -> POINT(-179.96 89.99) -> POINT(89.99 -179.96)
(3 0) -> POINT(-179.94 89.99) -> POINT(89.99 -179.94)
(4 0) -> POINT(-179.93 89.99) -> POINT(89.99 -179.93)
(5 0) -> POINT(-179.91 89.99) -> POINT(89.99 -179.91)
(6 0) -> POINT(-179.89 89.99) -> POINT(89.99 -179.89)
(7 0) -> POINT(-179.88 89.99) -> POINT(89.99 -179.88)
(8 0) -> POINT(-179.86 89.99) -> POINT(89.99 -179.86)
(9 0) -> POINT(-179.84 89.99) -> POINT(89.99 -179.84)
(10 0) -> POINT(-179.82 89.99) -> POINT(89.99 -179.82)