GeoJSON como objetos python y las librerías json y OGR

GeoJSON es un formato estándar abierto diseñado para representar elementos geográficos sencillos, junto con sus atributos no espaciales, basado en JavaScript Object Notation. El formato es ampliamente utilizado en aplicaciones de cartografía en entornos web al permitir el intercambio de datos de manera rápida, ligera y sencilla.

La librería OGR contiene el driver que permite manipular directamente este tipo de objetos desde el punto de vista geográfico. Sin embargo, se va a usar el módulo json para observar cómo es la sintaxis interna de estos archivos que pueden ser abiertos con un editor de texto convencional. Por ejemplo, el siguiente achivo geojson (Points.geojson), contiene una serie de 10 puntos con una proyección correspondiente a un EPSG:32612 (UTM12N/WGS84).

geojson1

Técnicamente, es una lista de diccionarios que, con el método ‘load’ del módulo json, podemos acceder como variables conociendo los keys correspondientes. Nuestro deseo es usar los métodos ‘CoordinateTransformation’ y ‘Transform’ de OGR para reproyectar los puntos desde EPSG:32612 a EPSG:4326 (long/lat WGS84). Estos métodos de OGR requieren una geometría GeoJSON y con el módulo json tenemos que aplicar el método ‘dumps’ para transformar la sintaxis de diccionario en un string WKT con la “geometría”.

Sin embargo, el string con la “geometría” es simplemente eso: un string. No es un objeto de ogr.Geometry. Para ello hay que emplear el método ‘CreateGeometryFromJson’ para obtenerla.

El script a continuación permite obtener las geometrías de los 10 features del archivo ‘GeoJSON’ considerado anteriormente, reproyectarlas desde 32612 a 4326, imprimir ambas geometrías en formato WKT y para la reproyectada sólo en formato QgsPoint. Finalmente, se crea una memory layer con los rasgos reproyectados.

from osgeo import ogr, osr
import os, json

# get the input layer
inDataSet = "/home/zeito/pyqgis_data/Points.geojson"

with open(inDataSet) as f:
    data = json.load(f)

try:
    srid = data['crs']['properties']['name'].split("::")[1]
    print srid

except KeyError:
    print "There is not crs assigned"

# input SpatialReference
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(int(srid))

# output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

geoms = []

for feature in data['features']:
    geom = feature['geometry']
    geom = json.dumps(geom)
    point = ogr.CreateGeometryFromJson(geom)
    print point
    # reproject the geometry
    point.Transform(coordTrans)
    print point
    print QgsPoint(point.GetX(), point.GetY())
    geoms.append(point.ExportToWkt())

epsg = 4326

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(geoms)) ]

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

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

La ejecución del código anterior imprime lo siguiente en la Python Console de QGIS:

32612
POINT (356348.441770963894669 4471883.137573418207467)
POINT (-112.692404143117344 40.385165903034)
(-112.692,40.3852)
POINT (355923.20246985956328 4472242.672360489144921)
POINT (-112.697493242224596 40.388330180833222)
(-112.697,40.3883)
POINT (356076.942377963103354 4472561.609190793707967)
POINT (-112.695754708141038 40.391228854998865)
(-112.696,40.3912)
POINT (356361.677787507011089 4473127.58941429015249)
POINT (-112.692528881768581 40.396374789806792)
(-112.693,40.3964)
POINT (354789.697379292920232 4472720.755257323384285)
POINT (-112.710951402468609 40.392438615370104)
(-112.711,40.3924)
POINT (355338.375289986783173 4473653.651575320400298)
POINT (-112.70470120130021 40.400934981984285)
(-112.705,40.4009)
POINT (355306.925854397937655 4472109.033362134359777)
POINT (-112.70472073826646 40.387019936976053)
(-112.705,40.387)
POINT (355931.519312235992402 4473130.85433157440275)
POINT (-112.697596186781013 40.396329900013164)
(-112.698,40.3963)
POINT (356272.675103192566894 4472649.188664698041975)
POINT (-112.693469229885366 40.392051321467449)
(-112.693,40.3921)
POINT (356816.558795929304324 4473194.782385258004069)
POINT (-112.687186211632834 40.397058200040007)
(-112.687,40.3971)

y permite visualizar los rasgos reproyectados en la Map View de QGIS como se expone a continaución:

geojson2

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

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