Generación de GeoJSON con Python

Question

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

Esto es fácil de hacer 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()

Como tengo todas las geometrías en un diccionario por geocódigo (self.geomdict) simplemente creo las características, establezco los campos y clono las geometrías de la capa preexistente (se omite el código de carga de esa capa para simplificar). 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. como desde el mapa de origen)

¿Cómo se exporta la colección de características generada como en el caso anterior?

Answer 1

Si tienes un entorno de desarrollo GDAL/OGR (cabeceras, librerías), podrías simplificar radicalmente tu código utilizando Fiona . Leer características de un shapefile, añadir nuevos atributos y escribirlos como GeoJSON es sólo 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)

Answer 2

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

fe.ExportToJson()

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

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

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

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

• Un CRS con nombre (por ejemplo, un URN 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 un dolor para obtener de OGR. En su lugar, si escribimos la proyección en un archivo en el disco que podemos referenciar con el URI. Podemos coger la proyección del objeto 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"
                                 }
                             }

Answer 3

Este es el más simple y fácil en Fiona. puede establecer el SRS para la salida GeoJSON.

import fiona
from fiona.crs import from_epsg

source= fiona.open('shp/second_shp.shp', 'r', encoding = 'utf-8')

with fiona.open('tool_shp_geojson/geojson_fiona.json','w',  driver ="GeoJSON", schema=source.schema, encoding = 'utf-8', crs=fiona.crs.from_epsg(4326)) as geojson:
     geojson.write(feat)