Shapefile de puntos en las celdas del perímetro de una máscara ráster

Supongamos que tenemos rasgos (features) de un archivo vectorial sobre celdas que albergan ‘nodata values’ o estén más allá de la extensión de un ráster y, aún así, queremos asociarlas a algún valor del ráster. Lo más lógico es que sea el que corresponda a la celda de éste más cercana al punto.

El problema no es muy complicado si el ráster no tiene ‘nodata values’ porque, para los puntos más allá de su extensión, se puede agenciar fácilmente un algoritmo usando los vértices superior izquierdo e inferior derecho del ráster. El problema es más complejo cuando se tiene una máscara del ráster porque, además de los puntos fuera de la extensión, pueden existir puntos asociados a ‘nodata values’ y cada situación se convierte en particular.

Para el ráster de la imagen inferior voy a extraer la máscara con base en el shapefile que allí se observa. Los ‘nodata values’ corresponderán a -999.

mask1

Elresultado es el siguiente:

mask2

donde la bounding box (cuadrado negro) corresponde al polígono de extensión del ráster y que permite visualizar el área útil del ráster.

El código que se presenta a continuación contiene el algoritmo que permite encontrar los índices de (fila, columna) de todas las celdas del perímetro del ráster. Estos índices se utilizan posteriormente para encontrar las coordenadas de los puntos que corresponden a la parte central de dichas celdas y conformar con ellos una memory layer de tales puntos.

from osgeo import gdal

dataset = gdal.Open('/home/zeito/pyqgis_data/aleatorio_mask.tif')

band = dataset.GetRasterBand(1)
noData = band.GetNoDataValue()

print noData

cols = band.XSize
rows = band.YSize

print rows, cols

transform = dataset.GetGeoTransform()

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

data = band.ReadAsArray(0, 0, cols, rows)

idxs = []

try:
    for i, row in enumerate(data):
        for j, element in enumerate(row):
            if data[i-1][j] == noData and element != noData or \
               data[i+1][j] == noData and element != noData or \
               data[i][j-1] == noData and element != noData or \
               element != noData and data[i][j+1] == noData:
               idxs.append([i, j])

except IndexError:
    print "Done!"

x0 = xOrigin + pixelWidth/2
y0 = yOrigin - pixelHeight/2

points = []

for item in idxs:
    x = x0 + item[1]*pixelWidth
    y = y0 - item[0]* pixelHeight
    points.append(QgsPoint(x,y))

wkt = dataset.GetProjection()

crs = QgsCoordinateReferenceSystem()
crs.createFromWkt(wkt)

epsg = crs.postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Ejecutado el código anterior en la Python Console de QGIS se obtiene:

mask3

Para probar que el código también funcionaba aceptablemente para rásteres de mayor envergadura, se empleó con uno de 502 filas por 546 columnas; tal como se evidencia a continuación:

mask4

A continuación se tiene un ‘Zoom In’ de una zona inferior del ráster donde se observa que el perímetro de puntos fue conformado de la manera esperada.

mask5

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

Una respuesta a Shapefile de puntos en las celdas del perímetro de una máscara ráster

  1. Pingback: Asignación de valores ráster más cercanos a puntos sobre ‘None’ values con 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