如何在两个具有不同坐标参考系的地理空间栅格数据之间找到对应的像素？

Question

我有几张 CRS = EPSG:32610 的 Landsat 8 图像。我有一个地面实况图像，其中每个像素代表该区域的 class，CRS = EPSG:4326。

我试图在光栅化中转换 ground truth 的 CRS，documentation 是这样说的：

dst_crs = raster.crs
transform, width, height = calculate_default_transform(
    labels.crs, dst_crs, labels.width, labels.height, *labels.bounds)
kwargs = labels.meta.copy()
kwargs.update({
    'crs': dst_crs,
    'transform': transform,
    'width': width,
    'height': height
})

with rasterio.open('groundtruth.tif', 'w', **kwargs) as dst:
    reproject(
        source=rasterio.band(labels, 1),
        destination=rasterio.band(dst, 1),
        src_transform=labels.transform,
        src_crs=labels.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.nearest)

这是转换前的标签（ground truth）元数据：

    {'driver': 'GTiff',
     'dtype': 'uint8',
     'nodata': None,
     'width': 4268,
     'height': 4782,
     'count': 1,
     'crs': CRS.from_epsg(4326),
     'transform': Affine(0.00032590510777881894, 0.0, -122.486874,
                         0.0, -0.0003259224173985775, 40.110855)}

这是改造后的样子：

    {'driver': 'GTiff',
     'dtype': 'uint8',
     'nodata': None,
     'width': 3721,
     'height': 5316,
     'count': 1,
     'crs': CRS.from_epsg(32610),
     'transform': Affine(32.8370217224448, 0.0, 543729.4708124083,
                         0.0, -32.8370217224448, 4441798.7147279335)}

这是 Landsat 图像元数据：

{'driver': 'GTiff',
 'dtype': 'int16',
 'nodata': None,
 'width': 7463,
 'height': 7702,
 'count': 6,
 'crs': CRS.from_epsg(32610),
 'transform': Affine(30.0, 0.0, 504045.0,
        0.0, -30.0, 4421925.0)}

我知道我可以像这样在地球表面的图像中找到像素的位置：

raster.transform * (x , y)

对于图像中的每个 (x, y) 位置。我的问题是 Landsat 图像中的像素比例为 30 米。当我进行转换时，比例不是 30x30。如果我使用 labels.transform * (x, y) 它永远不会匹配 Landsat 图像中 raster.transform * (x, y) 中的任何像素，并且数字有很多小数。

如何找到对应的像素点？图片和真实图片来自同一站点，但未完全对齐且大小不同。

这种方法对于找到对应点是否正确？我的专业不是遥感，我不熟悉CRS的专业术语和细节。如果您尽可能简单地向我解释，我将不胜感激。

Answer 1

在您的用例中，您希望 labels 图像与 raster 图像完美对齐。

为此，您不应使用 calculate_default_transform 函数，因为正如您所注意到的，它给出的 transform、width 和 height 并不完美匹配您的参考 raster 图片。

相反，您应该直接将 raster 的 transform、width 和 height 变量传递给 reproject 函数：

kwargs = labels.meta.copy()
kwargs.update({
    'crs': raster.crs,
    'transform': raster.transform,
    'width': raster.width,
    'height': raster.height
})

with rasterio.open('groundtruth.tif', 'w', **kwargs) as dst:
    reproject(
        source=rasterio.band(labels, 1),
        destination=rasterio.band(dst, 1),
        src_transform=labels.transform,
        src_crs=labels.crs,
        dst_transform=raster.transform,
        dst_crs=raster.crs,
        resampling=Resampling.nearest)

完成此操作后，您的 labels 和 raster 图像应该完全对齐，并且您应该能够将两个图像的像素一一对应。

如何在两个具有不同坐标参考系的地理空间栅格数据之间找到对应的像素？

How to find corresponding pixels between two geospatial raster data with different Coordinate Reference System?

python

coordinate-systems

coordinate-transformation

rasterio