0

I have a few Landsat 8 images that have CRS = EPSG:32610. I have a ground truth image in which each pixel represents a class for that region with CRS = EPSG:4326.

I tried to convert the CRS of ground truth in rasterio as documentation says like this:

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)

this is the labels (ground truth) meta data before transformation:

    {'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)}

and this is after transformation:

    {'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)}

this is the Landsat image meta data:

{'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)}

I know that I can find the location of a pixel in an image on the earth surface like this:

raster.transform * (x , y)

for every (x, y) position in the image. My problem is that the scale of pixels in Landsat images is 30 meters. When I do the transformation the scale is not quite 30x30. If I use the labels.transform * (x, y) it never matches any pixel in raster.transform * (x, y) in the Landsat image, and the numbers have a lot of decimals.

How can I find the corresponding pixels? Images and ground truth image are from the same site but not exactly aligned and have different sizes.

Is this approach even correct for finding corresponding points? My major is not remote sensing and I am not familiar with jargon and details of CRS. I appreciate it if you please explain to me as simple as you can.

guampi
  • 306
  • 1
  • 8
leo
  • 802
  • 6
  • 15

1 Answers1

2

In your use case you want your labels image to perfectly align with your raster image.

To do that, you should not use the calculate_default_transform function because as you have noticed, it gives a transform, width and height that do not perfectly match those of your reference raster image.

Instead, you should pass directly the transform, width and height variables of raster to the reproject function:

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)

Once you have done this, your labels and raster images should perfectly align, and you should be able to put in correspondence the pixels of both images one to one.

guampi
  • 306
  • 1
  • 8