0

What I want to do:

I have been trying to:

  • Open a raster
  • Make a copy
  • Rotate the copy by an angle with the raster center as a pivot
  • Keep the resolution constant in x/y dimension

Solution 1 with gdal:

#!/usr/bin/python
from optparse import OptionParser
import rasterio
from affine import Affine  # For easly manipulation of affine matrix
import scipy.ndimage
from rasterio.plot import reshape_as_raster, reshape_as_image
import numpy as np
from matplotlib import pyplot


def get_center(dataset):
    """This function return the pixel coordinates of the raster center
    """
    # We get the size (in pixels) of the raster using gdal
    #width, height = raster.RasterXSize, raster.RasterYSize
    width, height = dataset.width, dataset.height
    # We calculate the middle of raster
    xmed = width // 2
    ymed = height // 2
    return (xmed, ymed)


def rotate_geotransform(affine_matrix, angle, pivot):
    """This function generate a rotated affine matrix
    """
    affine_dst = affine_matrix * affine_matrix.rotation(angle, pivot)
    return(affine_dst)


def rotate(inputRaster, angle, outputRaster=None):
    outputRaster = 'rotated.TIF' if outputRaster is None else outputRaster

    src_dataset = rasterio.open(inputRaster)
    # this is a 3D numpy array, with dimensions [band, row, col]
    Z = src_dataset.read()

    # raster rotation
    old_affine_matrix = src_dataset.transform
    pivot = get_center(src_dataset)
    new_affine_matrix = rotate_geotransform(old_affine_matrix, angle, pivot)

    # array rotation
    rotated_Z = scipy.ndimage.rotate(Z, angle, order=1, reshape=True, axes=(1,2), cval=np.nan)
    print(Z.shape)
    pyplot.imshow(reshape_as_image(Z))
    pyplot.show()
    print(rotated_Z.shape)
    pyplot.imshow(reshape_as_image(rotated_Z))
    pyplot.show()

    new_dataset = rasterio.open(
        outputRaster,
        'w',
        driver='GTiff',
        height=rotated_Z.shape[1],
        width=rotated_Z.shape[2],
        count=rotated_Z.shape[0],
        dtype=Z.dtype,
        crs=src_dataset.crs,
        transform=new_affine_matrix
    )
    new_dataset.write(rotated_Z)
    new_dataset.close()


def main(argv):
    parser = OptionParser()
    parser.add_option("-o", "--output", type="str", dest="output", help="Rotated output raster name")
    (options, args) = parser.parse_args(argv)
    return rotate(args[0], float(args[1]), options.output)

if __name__ == '__main__':
    import sys
    main(sys.argv[1:])

This seems to rotate the first figure into the second: raster 1 rotated raster by 60 degree

Problem of solution 1

I am not sure my code does exactly what I intend it to do. Visually it seems ok, but is the coordinate reference system correctly calculated? It does not seem to keep the resolution constant (divided by appr. 2 in the rotated raster).

Solution 2 with rasterio and Affine packages:

#!/usr/bin/python
from optparse import OptionParser
import rasterio
from affine import Affine  # For easly manipulation of affine matrix
from rasterio.warp import reproject, Resampling
import numpy as np

def get_center(dataset):
    """This function return the pixel coordinates of the raster center
    """
    width, height = dataset.width, dataset.height
    # We calculate the middle of raster
    x_pixel_med = width // 2
    y_pixel_med = height // 2
    # The convention for the transform array as used by GDAL (T0) is to reference the pixel corner
    T0 = dataset.transform
    # We want to instead reference the pixel centre, so it needs to be translated by 50%:
    T1 = T0 * Affine.translation(0.5, 0.5)
    # to transform from pixel coordinates to world coordinates, multiply the coordinates with the matrix
    rc2xy = lambda r, c: T1 * (c, r)
    # get the coordinates for a raster in the first row, second column (index [0, 1]):
    return rc2xy(y_pixel_med, x_pixel_med)

def rotate(inputRaster, angle, outputRaster=None):
    outputRaster = 'rotated.tif' if outputRaster is None else outputRaster
    ### Read input
    source = rasterio.open(inputRaster)

    ### Rotate the affine
    pivot = get_center(source)
    pixel_size_x, pixel_size_y = source.res
    print("\nPivot coordinates:", pivot)
    new_transform = source.transform * Affine.rotation(angle, pivot) * Affine.scale(1)
    # this is a 3D numpy array, with dimensions [band, row, col]
    Z_source = source.read(masked=True)
    # Create destination raster
    destination = rasterio.open( outputRaster, 'w',
        driver='GTiff',
        height=source.height,
        width=source.width,
        count=source.count,
        crs=source.crs,
        dtype=Z_source.dtype,
        nodata=source.nodata,
        transform=new_transform)
    # Reproject pixels
    dst_shape = (destination.count, destination.height, destination.width)
    Z_destination = np.empty(dst_shape)
    Z_destination[:] = source.nodata

    reproject(
        Z_source,
        Z_destination,
        src_transform=source.transform,
        src_crs=source.crs,
        dst_transform=destination.transform,
        dst_crs=destination.crs,
        resampling=Resampling.average)

    destination.write(Z_destination)
    source.close()
    destination.close()
    return

def main(argv):
    parser = OptionParser()
    parser.add_option("-o", "--output", type="str", dest="output", help="Rotated output raster name")
    (options, args) = parser.parse_args(argv)
    return rotate(args[0], float(args[1]), options.output)

if __name__ == '__main__':
    import sys
    main(sys.argv[1:])

But somehow, the rotation seems to send the raster quite far, so I guess there is a bug with the code that should use the raster center cell as a pivot. However, I can't find what my mistake is.

3
  • 1
    I think the problem is related to the reshape parameter of the scipy.ndimage.rotate function. We can see in the documentation that this parameter changes the resolution of the image so that it keeps the same number of rows and columns. However I don't understand why you use the rotate function to perform the rotation, you can do it directly with rasterio or gdal.
    – Atm
    Jan 4, 2022 at 10:31
  • Thank you for you information about the reshape parameter. I did try to use rasterio and gdal (in the last case using the same resources you pointed to). But I did not succeed to visualize them rotated, and somehow that's how I ended up using the rotate function. If will try to use the rasterio issue you linked to and come back if I can't do any progress.
    – WaterFox
    Jan 11, 2022 at 16:23
  • @Atm I came with a code that uses rasterio (see Solution 2)! But my pivot seems off :( Any idea of what I got wrong?
    – WaterFox
    Jan 13, 2022 at 0:12

1 Answer 1

0

So, I got it working using rasterio. The problem was that the Affine pivot parameter is actually in pixel coordinates, and I was passing real world coordinates. So few lines of code later, this is rotating as exected!

#!/usr/bin/python
from optparse import OptionParser
import rasterio
from affine import Affine  # For easly manipulation of affine matrix
from rasterio.warp import reproject, Resampling
import numpy as np

def get_center_pixel(dataset):
    """This function return the pixel coordinates of the raster center
    """
    width, height = dataset.width, dataset.height
    # We calculate the middle of raster
    x_pixel_med = width // 2
    y_pixel_med = height // 2
    return (x_pixel_med, y_pixel_med)

def rotate(inputRaster, angle, scale=1, outputRaster=None):
    outputRaster = 'rotated.tif' if outputRaster is None else outputRaster

    ### Read input
    source = rasterio.open(inputRaster)
    assert source.crs == 'EPSG:4326', "Raster must have CRS=EPSG:4326, that is unprojected lon/lat (degree) relative to WGS84 datum"

    ### Rotate the affine about a pivot and rescale
    pivot = get_center_pixel(source)
    #pivot = None
    print("\nPivot coordinates:", source.transform * pivot)
    new_transform = source.transform * Affine.rotation(angle, pivot) * Affine.scale(scale)

    # this is a 3D numpy array, with dimensions [band, row, col]
    data = source.read(masked=True)
    kwargs = source.meta
    kwargs['transform'] = new_transform

    with rasterio.open(outputRaster, 'w', **kwargs) as dst:
        for i in range(1, source.count + 1):
            reproject(
                source=rasterio.band(source, i),
                destination=rasterio.band(dst, i),
                src_transform=source.transform,
                src_crs=source.crs,
                dst_transform=new_transform,
                dst_crs=dst.crs,
                resampling=Resampling.average)
    return

def main(argv):
    parser = OptionParser()
    parser.add_option("-o", "--output", type="str", dest="output", help="Rotated output raster name")
    (options, args) = parser.parse_args(argv)
    return rotate(args[0], float(args[1]), float(args[2]), options.output)

if __name__ == '__main__':
    import sys
    main(sys.argv[1:])

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.