# 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]

# 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:

# 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
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]
# 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.

• 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. 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? Jan 13, 2022 at 0:12

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

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]
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:])
``````