2

I am struggling with a task that seems like it should be quite simple. I have two rasters with the same size, resolution, CSR, but slightly different bounds. I need to perform a fairly simple calculation across, which requires that I align the two rasters. If I understand correctly, I need to perform some sort of reprojection on one of the rasters to align it with the other.

I see that there is function that provides a transform from the bounds, rasterio.transform.from_bounds(west, south, east, north, width, height), which seems like it would be the easiest way to go about this. Based on what I've found online, it seems like I need to do something like following to define and then apply the appropriate transform:

with rasterio.open('raster_1.tif') as src:
    west, south, east, north = src.bounds
    width = src.width'
    height = src.height
    transform = from_bounds(west, south, east, north, width, height)
    kwargs = src.meta.copy()
    kwargs.update({

        'crs': src.crs,
        'transform': transform,
        'width': width,
        'height': height})


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

When I check the bounds on the second raster, I see that it does indeed now line up with the first, except it seems that I have totally changed all of the data in raster. What have done wrong here?

3
  • 2
    If both rasters are already correctly georeferenced then it sounds like you just need to read the same geographic area from each of them. You could do this with a windowed read, creating a single window from bounds then reading that window from each raster
    – mikewatt
    Commented Jul 21, 2021 at 23:49
  • This looks promising, thank you! I'm getting the message that the window from bounds function requires a transform object to calculate the window. What sort of transform do I need if I don't actually want to warp the raster?
    – Max Miller
    Commented Jul 22, 2021 at 2:41
  • 2
    You need two windows (with same bounds), one per raster. You use each rasters dataset.transform as the transform for thw windpw from bounds function..
    – user2856
    Commented Jul 22, 2021 at 4:11

1 Answer 1

0

I have created align_extent_raster function to align the extent of one raster file (GeoTIFF) with another raster file (GeoTIFF), and then save the aligned raster to a specified output path.

The function parameters are:

  • raster_reference: The file path to the reference raster.
  • tiff_to_be_aligned: The file path to the raster that needs to be aligned to match the reference raster.
  • output_path: The file path where the aligned raster will be saved.
def align_extent_raster(raster_reference, tiff_to_be_aligned, output_path):
    with rasterio.open(raster_reference) as src1:
        data1 = src1.read(1)
        profile1 = src1.profile
        bounds1 = src1.bounds
        transform_1 = src1.transform
        no_data_1 = src1.nodata
        

    with rasterio.open(tiff_to_be_aligned) as src2:
        data2 = src2.read(1)
        profile2 = src2.profile
        bounds2 = src2.bounds
        transform_2 = src2.transform

    window = rasterio.windows.from_bounds(*bounds1, transform=transform_2)

    with rasterio.open(tiff_to_be_aligned) as src2:
        smaller_array = src2.read(window=window, boundless=True) 

    profile1.update(nodata= no_data_1)
    
    # Export it to another image 
    with rasterio.open(output_path, 'w', **profile1) as dst:
        dst.write(smaller_array)

Explanation

Both images are read and metadata are saved.

window = rasterio.windows.from_bounds(*bounds1, transform=transform_2)

The above-mentioned code code creates a window using the bounds of the reference raster (bounds1).

with rasterio.open(tiff_to_be_aligned) as src2:
    smaller_array = src2.read(window=window, boundless=True)

Then, the function reads the data from the raster to be aligned within the window that corresponds to the bounds of the reference raster and saves it to smaller_array. The last step is to save the aligned image.

Images

Reference image (red palette)

enter image description here

Reference image (red palette) overlaid by to be aligned image (blue palette) enter image description here

Reference image (red palette) overlaid by the aligned image (blue palette) after running the function

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.