sábado, 16 de noviembre de 2019

Cómo rotar un ráster con el módulo python de GDAL/OGR

En este hilo de gis.stackexchange.com se refiere la necesidad de rotar un ráster 90 grados en sentido antihorario debido a que los productos de 3B43 TRMM Multisatellite Precipitation monthly están supuestamente rotados 90 grados en sentido horario. Asumiendo que esto es cierto, con la ayuda de un algoritmo encontrado aquí, se desarrolló el codigo siguiente para rotar el array 90 grados en sentido antihorario y luego producir el raśter con los métodos de GDAL.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
from osgeo import gdal, osr
import numpy as np

dataset = gdal.Open("/home/zeito/pyqgis_data/test_0_5_0_5_clip.tif")
band = dataset.GetRasterBand(1)

array = dataset.ReadAsArray()

#rotate raster 90 degrees anticlockwise
array_rot = (list(list(x) for x in zip(*matrix))[::-1])

geotransform = dataset.GetGeoTransform()
wkt = dataset.GetProjection()

# Create gtif file
driver = gdal.GetDriverByName("GTiff")
output_file = "/home/zeito/pyqgis_data/rotated_raster.tif"

dst_ds = driver.Create(output_file,
                       band.XSize,
                       band.YSize,
                       1,
                       gdal.GDT_Float32)

array_rot = np.array(array_rot)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( array_rot )
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset

dataset = None
dst_ds = None

Aunque se probó con muchos píxeles, para ilustrar la comprobación se escogió el de la posición (16,19) en el ráster original. Este corresponde al valor de 803; tal como se desprende de la imagen siguiente con la ayuda del plugin Value Tool.


Este mismo valor, una vez realizada la rotación de 90 grados en sentido antihorario, debería estar en la posición (16,0) del ráster rotado. Es lo que se corrobora en la imagen siguiente:




No hay comentarios: