Generando GeoJSON con Python

14

Quiero crear un archivo GeoJSON mediante programación utilizando polígonos de un shapefile pero agregando atributos desde mi propia aplicación.

Esto se hace fácilmente para un shapefile:

def create_data_dayer(self,varlist, data):
    """
    Creates a new shape to contain data about nodes.
    varlist is the list of fields names associated with
    the nodes.
    data is a list of lists whose first element is the geocode
    and the remaining elements are values of the fields, in the
    same order as they appear in varlist.
    """
    if os.path.exists(os.path.join(self.outdir,'Data.shp')):
        os.remove(os.path.join(self.outdir,'Data.shp'))
        os.remove(os.path.join(self.outdir,'Data.shx'))
        os.remove(os.path.join(self.outdir,'Data.dbf'))
    # Creates a new shape file to hold the data
    if not self.datasource:
        dsd = self.driver.CreateDataSource(os.path.join(self.outdir,'Data.shp'))
        self.datasource = dsd
        dl = dsd.CreateLayer("sim_results",geom_type=ogr.wkbPolygon)
    #Create the fields
    fi1 = ogr.FieldDefn("geocode",field_type=ogr.OFTInteger)
    dl.CreateField(fi1)
    for v in varlist:
        #print "creating data fields"
        fi = ogr.FieldDefn(v,field_type=ogr.OFTString)
        fi.SetPrecision(12)
        dl.CreateField(fi)

    #Add the features (points)
    for n,l in enumerate(data):
        #Iterate over the lines of the data matrix.
        gc = l[0]
        try:
            geom = self.geomdict[gc]
            if geom.GetGeometryType() != 3: continue
            #print geom.GetGeometryCount()
            fe = ogr.Feature(dl.GetLayerDefn())
            fe.SetField('geocode',gc)
            for v,d in zip (varlist,l[1:]):
                #print v,d
                fe.SetField(v,str(d))
            #Add the geometry
            #print "cloning geometry"
            clone = geom.Clone()
            #print geom
            #print "setting geometry"
            fe.SetGeometry(clone)
            #print "creating geom"
            dl.CreateFeature(fe)
        except: #Geocode not in polygon dictionary
            pass
        dl.SyncToDisk()

ya que tengo todas las geometrías en un diccionario por geocodificación (self.geomdict) simplemente creo las entidades, configuro los campos y clono las geometrías de una capa preexistente (el código que carga esa capa se omite por simplicidad). Todo lo que necesito ahora es una forma de generar el GeoJSON a partir de la combinación de campos y geometrías, naturalmente con la ayuda de OGR para obtener el resto del archivo correcto (CRS, etc. a partir del mapa de origen)

¿Cómo exportar la colección de características generada como anteriormente?

    
pregunta fccoelho 19.11.2012 - 10:30

2 respuestas

13

Felizmente, OGR puede hacer esto por usted, ya que los objetos ogr.Feature y ogr.Geometry tienen los métodos ExportToJson() . En tu código;

fe.ExportToJson()

Y dado que los objetos geojson FeatureCollection son simplemente diccionarios con un type de FeatureCollection y un objeto features que contiene una lista de objetos Feature.

feature_collection = {"type": "FeatureCollection",
                      "features": []
                      }

feature_collection["features"].append(fe.ExportToJson())

El objeto CRS en una colección de características puede ser de dos tipos:

  • Un CRS denominado (por ejemplo, una URG de OGC o un código EPSG)
  • Un objeto de enlace con un URI y un tipo como "proj4"

Dependiendo de su formato de datos, es muy probable que el nombre sea difícil de obtener de OGR. En cambio, si escribimos la proyección en un archivo en el disco que podemos hacer referencia con la URI. Podemos tomar la proyección del objeto de capa (que tiene varias funciones de exportación)

spatial_reference = dl.GetSpatialRef()

with open("data.crs", "wb") as f:
    f.write(spatial_reference.ExportToProj4())

feature_collection["crs"] = {"type": "link",
                             "properties": {
                                 "href": "data.crs",
                                 "type": "proj4"
                                 }
                             }
    
respondido por el om_henners 19.11.2012 - 11:32
29

Si tienes un entorno GDAL / OGR dev (encabezados, libs), podrías simplificar radicalmente tu código usando Fiona . Para leer características de un shapefile, agregue nuevos atributos y escríbalos ya que GeoJSON es solo un puñado de líneas:

import fiona
import json

features = []
crs = None
with fiona.collection("docs/data/test_uk.shp", "r") as source:
    for feat in source:
        feat['properties'].update(...) # with your attributes
        features.append(feat)
    crs = " ".join("+%s=%s" % (k,v) for k,v in source.crs.items())

my_layer = {
    "type": "FeatureCollection",
    "features": features,
    "crs": {
        "type": "link", 
        "properties": {"href": "my_layer.crs", "type": "proj4"} }}

with open("my_layer.json", "w") as f:
    f.write(json.dumps(my_layer))
with open("my_layer.crs", "w") as f:
    f.write(crs)
    
respondido por el sgillies 20.11.2012 - 17:22

Lea otras preguntas en las etiquetas