Rásteres de longitud, latitud con GDAL/PyQGIS

En un post anterior se consideró la obtención de rásteres de longitud, latitud mediante la librería raster del lenguaje R. En esta ocasión se hará lo mismo pero con la librería GDAL y las clases de PyQGIS. Aunque el procedimiento de R parece exageradamente sencillo en comparación con este, hay que recordar que se puede hacer lo mismo encapsulando todo dentro de una función; sobre todo si el objetivo es crear un plugin.

Para probar el código se considerará un ráster relativamente extenso (1388 filas x 2105 columnas, 2.921.740 píxeles); tal como puede observarse en la imagen siguiente (pseudocolor, con rampa de colores ‘Spectral’, 5 clases):

longlat1

El código completo es el siguiente:

from osgeo import gdal, osr
import numpy as np
import os

layer = iface.activeLayer()
provider = layer.dataProvider()

my_path = provider.dataSourceUri()
(root,filename) = os.path.split(my_path)

dataset = gdal.Open(my_path)
band = dataset.GetRasterBand(1)

geotransform_src = dataset.GetGeoTransform()

extent = layer.extent()

xmin = extent.xMinimum()
ymax = extent.yMaximum()

rows = layer.height()
columns = layer.width()

xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

xinit = xmin + xsize/2
yinit = ymax - ysize/2

print "Processing..."

xcoor_values = []
ycoor_values = []

for i in range(rows):
    for j in range(columns):
        x = xinit + j*xsize 
        y = yinit
        xcoor_values.append(x)
        ycoor_values.append(y)
    xinit = xmin + xsize/2
    yinit -= ysize

driver = gdal.GetDriverByName("GTiff")

output_fileX = root + "/Xcoor_raster.tif"
output_fileY = root + "/Ycoor_raster.tif"

# Create GTif with rows and columns of base raster 
dst_dsx = driver.Create(output_fileX, 
                       columns, 
                       rows, 
                       1, 
                       gdal.GDT_Float32)  

dst_dsy = driver.Create(output_fileY, 
                       columns, 
                       rows, 
                       1, 
                       gdal.GDT_Float32) 

dst_srsx = osr.SpatialReference()
dst_srsy = osr.SpatialReference()

dst_srsx.ImportFromEPSG(4326)
dst_srsy.ImportFromEPSG(4326)

#Raster of x coordinates
raster_x = np.asarray(np.reshape(xcoor_values, (rows,columns)))

dst_dsx.GetRasterBand(1).WriteArray( raster_x )

dst_dsx.SetGeoTransform(geotransform_src)

dst_dsx.SetProjection(dst_srsx.ExportToWkt())

#Raster of y coordinates
raster_y = np.asarray(np.reshape(ycoor_values, (rows,columns)))

dst_dsy.GetRasterBand(1).WriteArray( raster_y )

dst_dsy.SetGeoTransform(geotransform_src)

dst_dsy.SetProjection(dst_srsy.ExportToWkt())

dataset = None
dst_dsx = None
dst_dsy = None

print "Done!"

y produce rásteres de longitud, latitud cuyos valores corresponden al centro de la celda del ráster base.

Cuando se ejecuta en la Python Console, los rásteres resultantes se obtienen en un tiempo relativamente breve y se cargan en el Map Canvas desde la ruta donde se crearon. En este caso lucen como a continuación (pseudocolor, con rampa de colores ‘Spectral’, 5 clases).

Ráster de longitudes:

longlat2

Ráster de latitudes:

longlat3

Esta entrada fue publicada en GDAL, PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Rásteres de longitud, latitud con GDAL/PyQGIS

  1. Pingback: Plugin para obtener rásteres de (longitud, latitud) con GDAL/PyQGIS | El Blog de José Guerrero

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