3

I've been trying to perform translation of a raster using the reproject function from rasterio.warp.

I've calculated the translation amount in X and Y from a series of "control points" placed manually on a shapefile. I'm using these values to update the affine transform of the source raster and keep everything else the same. Both input and output should be in UTM Zone 23S (EPSG:32723), same GSD, width and height. Here's my code:

import rasterio
from rasterio.transform import Affine
from rasterio.warp import reproject, Resampling
from pathlib import Path

input_raster_path = Path('path/to/raster.tif')
output_raster_path = input_raster_path.parent / (input_raster_path.stem + '_translated.tif')

with rasterio.open(input_raster_path, 'r') as src:

    dx, dy = -13.02168809006301, -2.439212051903208  # Values in meters, calculated earlier
    translation = Affine(1, 0, dx, 0, 1, dy)
    new_transform = translation * src.transform
    print(dx, dy)
    print(src.transform)
    print(new_transform)  # Making sure the new transform values make sense

    profile = src.profile.copy()
    profile.update(transform=new_transform)
    
    with rasterio.open(output_raster_path, 'w', **profile) as out:

        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(out, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=new_transform,
                dst_crs=src.crs,
                resampling=Resampling.nearest)

The above code results in a slight difference in the raster visually (probably due to rounding in the GSD and nearest resampling), but with no perceivable translation in relation to the original. I've also tried to multiply dx and dy by 10 to see if the values were just too small to notice, but the exact same output was created.

Note: The metadata on the new raster seems correct, I can inspect it on QGIS and get the correct extent and transform, but the pixels are not translated at all.

EDIT 1: I was unaware of the function Affine.translation but it produces the same translation matrix I've created directly by providing coefficients to the Affine constructor.

2
  • What's the raster's cell size?
    – mkennedy
    Commented Jan 7, 2021 at 4:38
  • Around 7 cm (0.07 m). I've exported a subsampled version directly from QGIS with a cell size of 0.4 m for quicker testing but the problem persists.
    – Guilherme
    Commented Jan 8, 2021 at 12:49

1 Answer 1

1

Returning after some time to share my findings. The reproject method cannot be used in this manner (or it can, but will not result in a raster translation as I was expecting). I've tested this again with a mockup raster. The following code should run as-is in an ipython notebook or as a standalone script:

This first part creates a mock single-band raster at a specific location with a cell size of 1 meter.

import rasterio
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt

img_data = np.uint8(np.random.rand(600, 800) * 255)
plt.imshow(img_data, cmap='gray')  # sanity check

meta = {
    'driver': 'GTiff', 
    'crs': 'epsg:32723',  # UTM CRS (coords in meters)
    'width': img_data.shape[1], 
    'height': img_data.shape[0], 
    'transform': rasterio.transform.from_origin(467415.32, 8882035.29, 1, 1),  # Sample transform with a cell size of 1x1 m
    'count': 1, 
    'dtype': img_data.dtype
    }
    
with rasterio.open('raster_a.tif', 'w', **meta) as dst:
    dst.write(img_data, 1)

The second part is very similar to my question's original code. It's responsible for applying the translation to the original transform and then simply reading and writing the band data to the new translated raster.

from rasterio.transform import Affine

input_raster_path = Path('raster_a.tif')
output_raster_path = input_raster_path.parent / (input_raster_path.stem + '_translated.tif')

with rasterio.open(input_raster_path, 'r') as src:

    dx, dy = -10, -10  # Values in meters, calculated earlier
    translation = Affine(1, 0, dx, 0, 1, dy)
    new_transform = translation * src.transform
    print(dx, dy)
    print(src.transform)
    print(new_transform)  # Making sure the new transform values make sense

    profile = src.profile.copy()
    profile.update(transform=new_transform)
    
    with rasterio.open(output_raster_path, 'w', **profile) as out:
        for i in range(1, src.count + 1):
            band_data = src.read(i)
            out.write(band_data, i)

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.