Reproyección punto a punto de un archivo vectorial con PyQGIS 3

En un post previo se consideró una forma de reproyectar un shapefile tipo polígono utilizando la versión 2 de PyQGIS. En éste, se va a modificar el código allí descrito para que pueda correr en PyQGIS 3 y, además, incluir todos los features del referido archivo vectorial.

En PyQGIS, las transformaciones entre diferentes sistemas de referencia espacial (CRS) se hacen a través de la clase QgsCoordinateTransform. El procedimiento sugiere crear los CRS fuente y destino y contruir la instancia de QgsCoordinateTransform con ellos. Entonces, se repite numerosas veces la llamada a la función xform() hasta culminar el proceso; tal como se presenta en el código a continuación. El shapefile usado presenta solo dos features (polígonos correspondientes a los estados Arizona y Utah de USA) a los cuales se les extraerá su geometría como polígono. Esta geometría viene dada como una lista de pares de coordenadas x,y para los anillos externos de cada polígono (polygon[0][0]). Por tanto, es sencillo implementar la función xform() para cada uno de estos puntos. Una vez obtenidos los puntos transformados entonces se recomponen en un nuevo shapefile tipo polígono.

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

layer = iface.activeLayer() #projecting_polygon.py polygons = [ feat.geometry().asMultiPolygon() for feat in layer.getFeatures() ] epsgSrc = layer.crs().postgisSrid() #source crs crsSrc = QgsCoordinateReferenceSystem(epsgSrc) #target crs crsDest = QgsCoordinateReferenceSystem(32612) xform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance()) n = len(polygons) polygons2 = [[] for i in range(n)] #anillo externo de cada poligono for i, polygon in enumerate(polygons): for point in polygon[0][0]: pt = xform.transform(point) polygons2[i].append(pt) polygons3 = [[] for i in range(n)] for i, item in enumerate(polygons2): polygons3[i].append(item) URI = "Polygon?crs=epsg:32612&field=id:integer""&field=area:double&index=yes" mem_layer = QgsVectorLayer(URI, "layer_memory", "memory") QgsProject.instance().addMapLayer(mem_layer) mem_layer.startEditing() feats = [ QgsFeature() for i in range(n) ] for i, feat in enumerate(feats): feat.setGeometry(QgsGeometry.fromPolygonXY(polygons3[i])) geom = feat.geometry() #Area determination area= geom.area() feat.setAttributes([i, area]) #set attributes values mem_layer.addFeature(feat) mem_layer.commitChanges() #stop editing and save changes

El código Python anterior refleja una transformación de WGS 84 en coordenadas geográficas a coordenadas UTM 12N (Utah parece que está comprendido totalmente en esa zona UTM). Una vez ejecutado y seleccionando el CRS con EPSG:32612, cuya proyección es en metros, los valores determinados para el atributo 'area' tienen sentido (219.740 km2 frente a los 219.887 km2 que reporta la Wikipedia; apenas un 0,07 % de diferencia).