Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Hi I want to programatically 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 ommited 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?

share|improve this question
add comment

2 Answers 2

up vote 6 down vote accepted

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"
                                 }
                             }
share|improve this answer
    
This is a good solution, because it doesn't add an extra dependency to my project like the (nice) solution of @sgillies –  fccoelho Nov 22 '12 at 15:42
    
I jus finished my testing with this solution and it worked nicely. However I had to handle manually when features had unicode characters in field names, since ogr.py didn't handle them properly. –  fccoelho Nov 24 '12 at 21:35
add comment

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)
share|improve this answer
1  
Fiona docs are killer! –  Chad Cooper Nov 20 '12 at 16:40
    
Thanks, Chad :) –  sgillies Nov 20 '12 at 20:22
    
Would vote up more than once if I could! –  om_henners Nov 20 '12 at 23:26
    
Isn't there a way to include the crs definition in GeoJSON? –  fccoelho Nov 22 '12 at 15:41
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.