Creando un ráster de distancia a la línea costera con PyQGIS 3 y GDAL/OGR

Con PyQGIS, cada centro de celda de un ráster puede ser calculado en un doble loop como un punto y, su distancia a la línea costera, evaluada mediante el método 'closestSegmentWithContext' de la clase QgsGeometry. Posteriormente, la librería GDAL puede ser usada para crear dicho ráster con estos valores de distancia en un array de filas x columnas.

El código siguiente fue empleado para crear un ráster de distancia (25 m × 25 m de resolución y 1000 filas x 1000 columnas) comenzando en el punto (397106.7689872353, 4499634.06675821); cerca de la costa oeste de USA (la proyección es UTM 10 N; EPSG:32610).

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 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

from osgeo import gdal, osr import numpy as np from math import sqrt registry = QgsProject.instance() line = registry.mapLayersByName('shoreline_10N') crs = line[0].crs() wkt = crs.toWkt() feats_line = [ feat for feat in line[0].getFeatures()] pt = QgsPoint(397106.7689872353, 4499634.06675821) xSize = 25 ySize = 25 rows = 1000 cols = 1000 raster = [ [] for i in range(cols) ] x = xSize/2 y = - ySize/2 for i in range(rows): for j in range(cols): point = QgsPointXY(pt.x() + x, pt.y() + y) tupla = feats_line[0].geometry().closestSegmentWithContext(point) raster[i].append(sqrt(tupla[0])) x += xSize x = xSize/2 y -= ySize data = np.array(raster) # Create gtif file driver = gdal.GetDriverByName("GTiff") output_file = "/home/zeito/pyqgis_data/distance_raster.tif" dst_ds = driver.Create(output_file, cols, rows, 1, gdal.GDT_Float32) #writting output raster dst_ds.GetRasterBand(1).WriteArray( data ) transform = (pt.x(), xSize, 0, pt.y(), 0, -ySize) #setting extension of output raster # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution dst_ds.SetGeoTransform(transform) # setting spatial reference of output raster srs = osr.SpatialReference() srs.ImportFromWkt(wkt) dst_ds.SetProjection( srs.ExportToWkt() ) dst_ds = None

Después de ejecutar el código anterior, empleando una capa vectorial de la línea costera, el ráster resultante fue cargado en QGIS y luce como en la imagen siguiente (pseudocolor con 5 classes y rampa Spectral).