Génération de GeoJSON avec Python

Je veux créer par programme un fichier GeoJSON en utilisant des polygones à partir d'un fichier de formes mais en ajoutant des attributs de ma propre application.

Cela se fait facilement pour un fichier de formes:

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

depuis que j'ai toutes les géométries sur un dictionnaire par géocodage (self.geomdict) je crée simplement les entités, définit les champs et clone les géométries à partir de la couche préexistante (code de chargement de cette couche omis pour plus de simplicité). Tout ce dont j'ai besoin maintenant est un moyen de générer le GeoJSON à partir de la combinaison de champs et de géométries, naturellement avec l'aide d'OGR pour obtenir le reste du fichier correctement (CRS, etc. à partir de la carte source)

Comment exporter la collection d'entités générée comme ci-dessus?

Heureusement, OGR peut le faire pour vous car les deux ogr.Featureet les ogr.Geometryobjets ont des ExportToJson()méthodes. Dans votre code;

fe.ExportToJson()

Et puisque GeoJSON FeatureCollection objets sont simplement des dictionnaires avec typede FeatureCollectionet un featuresobjet contenant une liste des objets de fonction.

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

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

L'objet CRS dans une collection d'entités peut être de deux types:

• Un CRS nommé (par exemple un URN OGC ou un code EPSG)
• Un objet lien avec un URI et un type tel que "proj4"

Selon votre format de données, il est très probable que le nom soit difficile à obtenir de OGR. Au lieu de cela, si nous écrivons la projection dans un fichier sur disque que nous pouvons référencer avec l'URI. Nous pouvons saisir la projection de l'objet calque (qui a plusieurs fonctions d'exportation)

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"
                                 }
                             }

Si vous avez un environnement de développement GDAL / OGR (en-têtes, bibliothèques), vous pouvez radicalement simplifier votre code en utilisant Fiona . Pour lire des entités à partir d'un fichier de formes, ajoutez de nouveaux attributs et écrivez-les car GeoJSON n'est qu'une poignée de lignes:

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)

C'est le plus simple et le plus facile de Fiona. vous pouvez définir le SRS pour la sortie 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)