El blog de José Guerrero: Algoritmo de Ramer–Douglas

El algoritmo de Ramer–Douglas–Peucker (RDP)¹, ² permite reducir el número de puntos utilizados en la aproximación de una curva. Se fundamenta en la búsqueda de puntos críticos a partir de una tolerancia lineal. La primera línea base corresponde al primer (anclaje) y último punto (flotante) de la línea original. Luego se calculan las distancias perpendiculares de todos los puntos intermedios. Si ninguna de estas distancias es mayor a la tolerancia sólo quedarán retenidos los puntos inicial y final de la línea. Sí se supera la tolerancia el punto con mayor distancia será retenido como punto crítico subdividiendo la línea original en dos secciones. En cada una se repite el proceso de manera independiente y así, sucesivamente, hasta que no haga falta efectuar subdivisión alguna de las líneas.

Una implementación del algoritmo en Python se encuentra ya en el módulo 'rdp' que puede disponibilizarse para la Python Console de QGIS a través de pip install. Una vez instalado, el módulo se probó con el código siguiente:

import processing
from rdp import rdp

layer = processing.getObjectFromName('argentina')
feats = [ feat for feat in layer.getFeatures() ]

for feat in feats:
    if feat.geometry().isGeosValid():
        print "geometry is valid"
    else:
        print "geometry is not valid for ID: ", feat.attribute("ID")
 
argentine = feats[0].geometry().asPolygon()[0]

print len(argentine)

res = rdp(argentine, epsilon = 0.1)

print len(res)

new_res = [ QgsPoint(point[0], point[1]) for point in res ]
new_geom_l = QgsGeometry.fromPolyline(new_res).exportToWkt()
new_geom_p = QgsGeometry.fromPolygon([new_res]).exportToWkt()

epsg = layer.crs().postgisSrid()

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

mem_layer_l = QgsVectorLayer(uri_l,
                             'mem_layer_line',
                             'memory')

prov = mem_layer_l.dataProvider()

feat = QgsFeature()

feat.setAttributes([0])
feat.setGeometry(QgsGeometry.fromWkt(new_geom_l))

prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer_l)

path = processing.runalg('saga:convertlinestopolygons',
                         'mem_layer_line',
                         None)

layer_polygon = QgsVectorLayer(path['POLYGONS'],
                               'layer_polygon',
                               'ogr')

QgsMapLayerRegistry.instance().addMapLayer(layer_polygon)

Utiliza el feature correspondiente al territorio continental de Argentina que fue obtenido del shapefile world_borders; disponible en Internet. Este contiene digitalizados 3450 puntos y como sé que el referido shapefile contiene algunas geometrías inválidas se hace una validación previa antes de la ejecución del algoritmo. La tolerancia empleada fue de 0.1 grados lo que simplifica el shapefile a sólo 210 puntos. Los resultados obtenidos se visualizan, en primer lugar, mediante una memory layer que sirve como entrada para obtener un capa tipo polígono mediante el método saga:convertlinestopolygons de Processing.

Después de ejecutar el algoritmo en la Python Console de QGIS se obtiene algo similar a lo siguiente:

El Zoom Full prácticamente hace indistinguibles los detalles. No obstante, un acercamiento con una menor escala permite observar claramente los efectos del algoritmo de simplificación en la imagen siguiente: