Cómo “cortar” un ráster con GDAL/Python pero sin gdal_warp o gdal_translate

La utilidad que usa por defecto QGIS para cortar un ráster es gdal_translate. No obstante, es posible utilizar el método ‘ReadAsArray’ de la librería GDAL de Python para seleccionar el bloque completo que nos interesa “cortar”. Sólo es necesario conocer las coordenadas de un punto (p1) en la celda del vértice superior izquierdo y de un punto (p2) en la celda del vértice inferior derecho del referido bloque para poder determinar los índices (fila, columna) respectivos.

Los índices de (columna [i], fila [j]) se determinan por estas formulas:

i1 = int((p1[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - p1[1] ) / pixelHeight)
i2 = int((p2[0] - xOrigin) / pixelWidth)
j2 = int((yOrigin - p2[1]) / pixelHeight)

lo que da lugar al número de columnas y filas del ráster resultante siguiente:

new_cols = i2-i1+1
new_rows = j2-j1+1

La selección de un bloque cuya primera celda es diferente a la (0,0) involucra también la determinación de un nuevo parámetro de geotransform; tal como se indica a continuación:

new_x = xOrigin + i1*pixelWidth
new_y = yOrigin - j1*pixelHeight
new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])

Para probar mi procedimiento con esta situación:

cut_raster1

y el código siguiente:

from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
filename = "c:/pyqgis_data/aleatorio.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

print cols, rows

transform = dataset.GetGeoTransform()

print transform

p1 = (355217.199739, 4473171.2377)
p2 = (355911.113396, 4472582.9196)

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

print xOrigin, yOrigin

i1 = int((p1[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - p1[1] ) / pixelHeight)
i2 = int((p2[0] - xOrigin) / pixelWidth)
j2 = int((yOrigin - p2[1]) / pixelHeight)

print i1, j1
print i2, j2

new_cols = i2-i1+1
new_rows = j2-j1+1

data = band.ReadAsArray(i1, j1, new_cols, new_rows)

print data

new_x = xOrigin + i1*pixelWidth
new_y = yOrigin - j1*pixelHeight

print new_x, new_y

new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])

# Create gtif file 
driver = gdal.GetDriverByName("GTiff")

output_file = "c:/pyqgis_data/cut_raster.tif"

dst_ds = driver.Create(output_file, 
                       new_cols, 
                       new_rows, 
                       1, 
                       gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( data )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(new_transform)

wkt = dataset.GetProjection()

# 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

su ejecución en la Python Console de QGIS da lugar a:

cut_raster2

Para ejecutar el código anterior, pero esta vez con el mismo ráster reproyectado en lat/long WGS 84 (EPSG: 4326), fue necesario realizar las modificaciones siguientes:

.
.
.
filename = "/home/zeito/pyqgis_data/aleatorio_4326.tif"
.
.
.
p1 = (-112.706287864, 40.3963533407)
p2 = (-112.70234697, 40.3932517699)
.
.
.
output_file = "/home/zeito/pyqgis_data/cut_raster2.tif"
.
.
.

cuyo resultado se visualiza a continuación:

cut_raster3

El corte, en ambos casos, se produjo con la alineación esperada.

Esta entrada fue publicada en Código Python, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Cómo “cortar” un ráster con GDAL/Python pero sin gdal_warp o gdal_translate

  1. pedro ls dijo:

    Hola, no se si es el modo adecuado para hacerte una consulta, pero si asi fuera, te agradezco de antemano. El asunto es que (en QGis) al exportar el shp que esta categorizado a kml no me respeta la categorizacion realizada ni lo colores, ni tamaños, tal vez no es posible hacer lo que quiero, siendo asi, existe alguna forma de crear mapas categorizados y que pueda exportarlo a kml. Gracias Pedro

    Date: Fri, 24 Jun 2016 09:42:30 +0000 To: ls_pedro@hotmail.com

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s