将笛卡尔像素数据数组转换为纬/经像素数据数组。

Question

将笛卡尔像素数据数组转换为纬/经像素数据数组。

3

我有一张图片（基本上，我获得了1024x1024像素的原始图像数据）和该图像中心像素的经纬度位置。

每个像素代表相同的固定像素比例尺寸（例如每个像素30米）。

现在，我想将图像绘制到使用坐标参考系统“EPSG：4326”（WGS84）的地图上。

当我仅定义图像的经纬度角落并根据“像素大小*像素比例尺”的计算以及将距离从中心点转换为每个角落的经纬度坐标来绘制它时，我认为，图像未正确绘制到地图上。通过术语“未正确绘制”，我的意思是，图像似乎被移动了，而且图像的内容不在我期望它们在地图位置上。

我认为这是因为我“混合”了像素比例尺和“EPSG：4326”坐标参考系统。

现在，根据我所给出的信息，我是否可以使用Geotools将整个像素矩阵从固定像素比例基础转换为“EPSG：4326”坐标参考系统中的新像素矩阵？当然，转换必须依赖于我已经给出的经纬度中心位置和像素比例尺。

我想知道是否使用类似于这样的东西可以指引我正确的方向：

MathTransform transform = CRS.findMathTransform(DefaultGeocentricCRS.CARTESIAN, DefaultGeographicCRS.WGS84, true);

  DirectPosition2D srcDirectPosition2D = new DirectPosition2D(DefaultGeocentricCRS.CARTESIAN, degreeLat.getDegree(), degreeLon.getDegree());
  DirectPosition2D destDirectPosition2D = new DirectPosition2D();
  transform.transform(srcDirectPosition2D, destDirectPosition2D);

  double transX = destDirectPosition2D.x;
  double transY = destDirectPosition2D.y;

  int kmPerPixel = mapImage.getWidth / 1024; // It is known to me that my map is 1024x1024km ...

  double x = zeroPointX + ((transX * 0.001) * kmPerPixel);
  double y = zeroPointY + (((transX * -1) * 0.001) * kmPerPixel);

我从另一个SO线程中得到了这段代码，并已经对其进行了一些修改，但仍然想知道这是否是解决我的问题的正确起点。

我只假设我的原始图像坐标参考系统属于DefaultGeocentricCRS.CARTESIAN类型。有人能确认一下吗？

接下来，使用Geotools来解决这种问题是否是正确的起点，还是我走错了路？

此外，我想补充说明这将用于一个相对安静的动态系统。因此，我的图像更新速率约为10Hz，需要相应地执行变换。再次问一下，我的初步想法是否导致了解决方案，还是您有其他解决方案来解决我的问题？

非常感谢， Kiamur

- Kiamur

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ian Turton · Accepted Answer

这并不像听起来那么简单。你实际上是在试图使用一个平面正方形来定义一个球体（严格来说是椭球体）上的区域。因此，没有一种“正确”的方法可以做到完美无缺，总会有些扭曲变形。如果不知道你的图片来自哪里，就无法给出确切的答案，但以下代码提供了三个可能的答案：

第一个和第二个使用GeoTools的GeodeticCalculator来计算角点，使用方位角和距离。这些是上面的蓝色“正方形”和绿色“正方形”。蓝色正在直接计算角落，而绿色正在计算边缘，并从交点推断出角落（这就是为什么它更加正方形）。

final int width = 1024, height = 1024;
GeometryFactory gf = new GeometryFactory();
Point centre = gf.createPoint(new Coordinate(0,51));
WKTWriter writer = new WKTWriter();
//direct method
GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
calc.setStartingGeographicPoint(centre.getX(), centre.getY());
double height2 = height/2.0;
double width2 = width/2.0;
double dist = Math.sqrt(height2*height2+width2 *width2);
double bearing = 45.0;
Coordinate[] corners = new Coordinate[5];
for (int i=0;i<4;i++) {
  calc.setDirection(bearing, dist*1000.0 );
  Point2D corner = calc.getDestinationGeographicPoint();
  corners[i] = new Coordinate(corner.getX(),corner.getY());
  bearing+=90.0;
}
corners[4] = corners[0];
Polygon bbox = gf.createPolygon(corners);
System.out.println(writer.write(bbox));

double[] edges = new double[4];
bearing = 0;
for(int i=0;i<4;i++) {
  calc.setDirection(bearing, height2*1000.0 );
  Point2D corner = calc.getDestinationGeographicPoint();

  if(i%2 ==0) {
    edges[i] = corner.getY();
  }else {
    edges[i] = corner.getX();
  }
  bearing+=90.0;
}

corners[0] = new Coordinate( edges[1],edges[0]);
corners[1] = new Coordinate( edges[1],edges[2]);
corners[2] = new Coordinate( edges[3],edges[2]);
corners[3] = new Coordinate( edges[3],edges[0]);

corners[4] = corners[0];
bbox = gf.createPolygon(corners);
System.out.println(writer.write(bbox));

另一种方法是将中心点转换为一个“更平”的投影，并使用简单的加法计算角落，然后反转变换。为此，我们可以使用OGC WMS规范定义的AUTO投影生成以我们的点为中心的正射投影，这给出了非常类似于蓝色方块的红色“方块”。

String code = "AUTO:42003," + centre.getX() + "," + centre.getY();
// System.out.println(code);
CoordinateReferenceSystem auto = CRS.decode(code);
// System.out.println(auto);
MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84,
    auto);
MathTransform rtransform = CRS.findMathTransform(auto,DefaultGeographicCRS.WGS84);
Point g = (Point)JTS.transform(centre, transform);

width2 *=1000.0;
height2 *= 1000.0;
corners[0] = new Coordinate(g.getX()-width2,g.getY()-height2);
corners[1] = new Coordinate(g.getX()+width2,g.getY()-height2);
corners[2] = new Coordinate(g.getX()+width2,g.getY()+height2);
corners[3] = new Coordinate(g.getX()-width2,g.getY()+height2);
corners[4] = corners[0];
bbox = gf.createPolygon(corners);
bbox = (Polygon)JTS.transform(bbox, rtransform);
System.out.println(writer.write(bbox));

使用哪种解决方案是个人品味的问题，取决于你的图像来源，但我认为红色或蓝色可能是最好的选择。如果你需要以10Hz的速度进行操作，则需要测试它们的速度，但我怀疑转换图像将成为瓶颈。

一旦你满意了边界框的设置，你可以使用以下方法将未引用的图像转换为地理参考覆盖范围：

GridCoverageFactory factory = CoverageFactoryFinder.getGridCoverageFactory(null);
GridCoverage2D gc = factory.create("name", image, new ReferencedEnvelope(bbox.getEnvelopeInternal(),DefaultGeographicCRS.WGS84));
String fileName = "myImage.tif";
AbstractGridFormat format = GridFormatFinder.findFormat(fileName);
File out = new File(fileName);
GridCoverageWriter writer = format.getWriter(out);
try {
  writer.write(gc, null);
  writer.dispose();
} catch (IllegalArgumentException | IOException e) {
  // TODO Auto-generated catch block
  e.printStackTrace();
}