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:

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.