Как извлечь значения широты и долготы из изображения GeoTIFF, используя Java? - PullRequest
Новая
       63

Как извлечь значения широты и долготы из изображения GeoTIFF, используя Java?

0 голосов
/ 01 июня 2018

У меня есть спутниковое изображение, на котором я должен работать.Я прочитал файл, используя ImageIO и JAI.Теперь я хочу найти значения широты и долготы, соответствующие каждому пикселю, который хранится в файле TIFF.Заранее спасибо.

1 Ответ

0 голосов
/ 01 июня 2018

Если у вас есть 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)
...