Cómo reproyectar un shapefile usando PyQGIS

Según el PyQGIS developer cookbook v. 2.0, 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 transform() hasta culminar el proceso.

Para probar esto se va a usar un shapefile con solo dos features (polígonos correspondientes a los estados Arizona y Utah de USA) de los cuales, a través de la Python Console, se seleccionará uno de ellos por id (método select para id=1; estado de Utah) y se le extraerá su geometría como polígono. Esta geometría en realidad viene dada como una lista de pares de coordenadas x,y (tal como aclara la documentación); a pesar de su apariencia de lista de tuplas cuando se imprime. Por tanto, es sencillo implementar la función transform() para cada uno de estos puntos. Una vez obtenidos los puntos transformados entonces se recomponen en un nuevo shapefile tipo polígono usando el mismo procedimiento de un post anterior.

El código Python 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) y se indica a continuación:

au=iface.activeLayer()

au.select(1)

feature=au.selectedFeatures()

geom=feature[0].geometry().asPolygon()

print "\noriginal points\n", geom

crsSrc = QgsCoordinateReferenceSystem(4326)
crsDest = QgsCoordinateReferenceSystem(32612)

xform = QgsCoordinateTransform(crsSrc, crsDest)

geom2=geom

n = len(geom2[0])

for i in range(n):
    pt = xform.transform(QgsPoint(geom[0][i]))
    geom2[0][i]=pt

print "\ntransformed points\n", geom2

mem_layer = QgsVectorLayer("Polygon?crs=epsg:32612&field=id:integer"
                                            "&field=area:double&index=yes",
                                            "layer_memory",
                                            "memory")

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

mem_layer.startEditing()

feature = QgsFeature()

feature.setGeometry(QgsGeometry.fromPolygon(geom2))

#Area determination 
geom = feature.geometry()
 
area= geom.area()
 
#set attributes values
feature.setAttributes([1, area])
 
mem_layer.addFeature(feature, True)
 
#stop editing and save changes
mem_layer.commitChanges()

y acá lo tenemos en el Editor de la Python Console donde se imprimen algunos resultados intermedios con propósito de verificación:

projection2

La imagen de partida en el Map Canvas es ésta:

projection4

y la resultante es esta otra:

projection3

Observen que como ahora la 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).

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

3 respuestas a Cómo reproyectar un shapefile usando PyQGIS

  1. Adhemar Conde dijo:

    Hola José,
    Felicidad por el trabajo que vienes haciendo. Desde que encontré tu blog me he interesado muchísimo más en las herramientas SIG libres, poco a poco voy utilizando herramientas como Grass Gis, R y GDAL.

    Hoy te escribo porque quiero pedirte ayuda con un problema que tengo.
    Estas últimas semanas he intentando realizar el cálculo de la reflectancia aparente en imágenes Landsat 5, con el paquete Landsat para R y con GrassGis.
    Ocurre que en teoría los valores deberían estar en el rango de 0 a 1, pero yo obtengo valores de -0 a 1, no se si esto es normal (he realizado el procedimiento en R y Grass), pero con estos valores tengo problemas al calcular el NDVI.
    En resumen más o menos ese es el problema que tengo. También he visto en otros comentarios en tu mismo blog, que otro compañero ha tenido más o menos el mismo problema y la verdad no sé si a podido solucionarlo.

    Te dejo la dirección de mi correo: adhemarfabio@gmail.com. Así podría explicarte con más detalle este problema, tal vez podría enviarte la linea de código que tengo para R, los metadatos de la imagen, etc., tal vez estoy cometiendo errores. Por supuesto, si es que tu aceptarías ayudarme. Estoy absolutamente seguro, que a ti no te va ha tomar mucho tiempo hallar el error.

    También, te pido disculpas si esta no es la forma de entrar en contacto contigo. Pero no he encontrado otra manera.

    Te envío muchos saludos desde Bolivia – La Paz.
    Y me gustaría muchísimo recibir una respuesta positiva de tu parte.

    Adhemar.

  2. Pingback: Objetos de la Clase QgsCoordinateReferenceSystem en 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