Generating GeoJSON with Python

Question

I want to programmatically create a GeoJSON file using polygons from a shapefile but adding attributes from my own application.

This is easily done for a 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()

since I have all the geometries on a dictionary by geocode (self.geomdict) I simply create the features, set the fields and clone the geometries from pre-existing layer (code loading that layer omitted for simplicity). All I need now is a way to generate the GeoJSON from the combination of fields and geometries, naturally with the help of OGR to get the rest of the file right (CRS, etc. as from the source map)

How do export the feature collection generated as above?

Answer 1

If you've got an GDAL/OGR dev environment (headers, libs), you could radically simplify your code by using Fiona. To read features from a shapefile, add new attributes, and write them out as GeoJSON is just a handful of lines:

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

Happily OGR can do this for you as both ogr.Feature and ogr.Geometry objects have ExportToJson() methods. In your code;

fe.ExportToJson()

And since geojson FeatureCollection objects are simply dictionaries with a type of FeatureCollection and a features object containing a list of Feature objects.

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

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

The CRS object in a feature collection can be one of two types:

• A named CRS (e.g. an OGC URN or an EPSG code)
• A link object with a URI and a type such as "proj4"

Depending on your data format it's quite likely that the name is going to be a pain to get from OGR. Instead if we write the projection to a file on disk which we can reference with the URI. We can grab the projection from the layer object (which has several Export functions)

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

This is the simplest and easiest one in Fiona. you can set the SRS for output 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)