How to create line/polygon shapefiles from geojson using gdal in Python?

Question

I have put this together after reading How to create polygon shapefile from a list of coordinates using python gdal/ogr?, and How to write Shapely geometries to shapefiles?

I am attempting to read a GeoJSON response and put together 3 different shapefiles; one for each geometric type that may be contained in the response. Since I'm not very good with gdal/ogr, I figured I'd include a lengthy chunk of the script for evaluation.

Methods I'm using to make lines and polygons...

def create_line(coords):
    line = ogr.Geometry(ogr.wkbLineString)
    for coord in coords:
        line.AddPoint(float(coord[0], float(coord[1]))
    return line.ExportToWkt()

def create_polygon(coords):
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for coord in coords:
        ring.AddPoint(float(coord[0], float(coord[1]))

    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    return poly.ExportToWkt()

Preparation and Loops for reading the geojson and building the features...

#prep files
outdir = "C:\\Temp"
p_shp = os.path.join(outdir, "point.shp")
l_shp = os.path.join(outdir, "line.shp")
a_shp = os.path.join(outdir, "area.shp")

driver = ogr.GetDriverByName("ESRI Shapefile")
point_data_source = driver.CreateDataSource(os.path.join(outdir, p_shp)
line_data_source = driver.CreateDataSource(os.path.join(outdir, l_shp)
area_data_source = driver.CreateDataSource(os.path.join(outdir, a_shp)

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# Point shapefile setup
p_layer = point_data_source.CreateLayer("point_layer", srs, ogr.wkbPoint)
field_testfield = ogr.FieldDefn("fld_a", ogr.OFTString)
field_testfield.SetWidth(50)
p_layer.CreateField(field_testfield)

# Line shapefile setup
l_layer = line_data_source.CreateLayer("line_layer", srs, ogr.wkbLineString)
field_testfield = ogr.FieldDefn("fld_a", ogr.OFTString)
field_testfield.SetWidth(50)
l_layer.CreateField(field_testfield)

# Area shapefile setup
a_layer = area_data_source.CreateLayer("area_layer", srs, ogr.wkbPolygon)
field_testfield = ogr.FieldDefn("fld_a", ogr.OFTString)
field_testfield.SetWidth(50)
a_layer.CreateField(field_testfield)    

# Do the json reading
for i in jsonObj['features']:
    geo = i.get('geometry')
    geo_type = geo.get('type')
    if geo_type == 'Point':  # If it's a point feature, do this
        [xcoord, ycoord, field_value] = ['', '', 'test']
        coords = geo.get('coordinates')
        xcoord, ycoord = float(coords[0]), float(coords[1])
        # Get properties and build field values here
        if xcoord and ycoord and field_value:
            feature = ogr.Feature(p_layer.GetLayerDefn())
            feature.SetField('fld_a', field_value)

            wkt = "POINT ({} {})".format(xcoord, ycoord)
            point = ogr.CreateGeometryFromWkt(wkt)
            feature.SetGeometry(point)
            p_layer.CreateFeature(feature)
            feature = None
            [xcoord, ycoord, field_value] = ['', '', '']

    elif geo_type == 'LineString':  # Or if it's a line feature, do this
        [line, field_value] = ['', 'test']
        linecoords = geo.get('coordinates')
        line = create_line(linecoords)
        if line and field_value:
            feature = ogr.Feature(l_layer.GetLayerDefn())
            feature.SetField('fld_a', field_value)

            wkt = "LINESTRING ({})".format(line)
            line_feat = ogr.CreateGeometryFromWkt(wkt)
            feature.SetGeometry(line_feat)
            l_layer.CreateFeature(feature)
            feature = None
            [line, field_value] = ['', '']

    elif geo_type == 'Polygon':  # Or if it's a polygon, do this
        [poly, field_value] = ['', 'test']
        polycoords = geo.get('coordinates')
        poly = create_polygon(polycoords)
        if poly and field_value:
            feature = ogr.Feature(a_layer.GetLayerDefn())
            feature.SetField('fld_a', field_value)

            wkt = "POLYGON ({})".format(poly)
            area = ogr.CreateGeometryFromWkt(wkt)
            feature.SetGeometry(area)
            a_layer.CreateFeature(feature)
            feature = None
            [poly, field_value] = ['', '']

    else:
        print('Could not discern geometry')

The Result (UPDATED)

I'm now using Fiona to do this in the way @gene suggested below:

for i in jsonObj['features']
geo = i.get('geometry')
geo_type = geo.get('type')
props = i.get('properties')
value_a = props.get('value_a', 'UNK')
value_b = props.get('value_b', 'UNK')
if geo_type == 'Point':
    schema = {'geometry': 'Point', 'properties': {
        'field_a': 'str:50',
        'field_b': 'str:50'}}
    with fiona.open(point_shapefile, 'w', 'ESRI Shapefile', schema) as layer:
        try:
            layer.write({'geometry': json.loads(geo), 'properties': {
                    'field_a': value_a,
                    'field_b': value_b}})
        except TypeError as te:
            print('Problem creating point feature, {}'.format(str(te)))

This yields:

Problem creating point feature, expected string or buffer

I am able to access and print all values out of the json response object jsonObj. However, the Fiona layer writing is failing.

why not use the gdal/ogr geojson driver to read the json in? — Ian Turton, Jul 06 '17 at 14:36
@iant that sounds nice; can you direct me to documentation for doing this in python? Googling didn't turn up anything on my end and the gdal python cookbook does not seem to contain anything like "Save GeoJSON to Output Layer." — auslander, Jul 06 '17 at 14:59
See https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#get-geometry-from-each-feature-in-a-layer but use GeoJSON instead of ESRI Shapefile in the driver = ogr.GetDriverByName("ESRI Shapefile") line — Ian Turton, Jul 06 '17 at 15:20
When I tried dataSource = driver.Open(jsonObj, 0) and then did for feature in layer: geom = feature.GetGeometryRef() \ print geom it yielded RuntimeError: not a string — auslander, Jul 06 '17 at 15:32
@iant nope that doesn't work either, even if i try to coerce it to str. Still get the not-a-string error. Even if I do for feature in layer: print("{}".format(str(feature))) it will give me a not a string. Somehow I think it doesn't want to open up r.json() (the json response from Requests). Maybe just tweaking what i have is the best option? — auslander, Jul 06 '17 at 15:54
What is the result of print {'geometry': json.loads(geo), 'properties': {'field_a': value_a,'field_b': value_b}}) ? — gene, Jul 10 '17 at 16:33
I isolated json.loads(geo) and that's where the TypeError is coming from. When I print geo, it prints it successfully as {u'type': u'Point', u'coordinates': [-90.24565, 38.63979]} Perhaps it is already being read correctly as a python dict? — auslander, Jul 10 '17 at 16:44
The result should be {'geometry': {u'type': u'Point', u'coordinates': [-90.24565, 38.63979]}, 'properties': {'field_b': 'UNK', 'field_a': 'UNK'}} — gene, Jul 10 '17 at 17:32

score 4 · Accepted Answer · answered Jul 07 '17 at 14:29

For creating an OGR geometry from GeoJSON, see Python GDAL/OGR Cookbook: Create Geometry from GeoJSON

point = """{"type":"Point","coordinates":[108420.33,753808.59]}"""
line = """{ "type": "LineString", "coordinates": [ [100.0, 0.0], [101.0, 1.0] ]}"""
poly = """{ "type": "Polygon","coordinates": [[ [100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0] ],[ [100.2, 0.2], [100.8, 0.2], [100.8, 0.8], [100.2, 0.8], [100.2, 0.2] ]]}"""
from osgeo import ogr
pointogr = ogr.CreateGeometryFromJson(point)
lineogr = ogr.CreateGeometryFromJson(line)
polyogr = ogr.CreateGeometryFromJson(poly)

After that

driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource('line.shp')
layer = ds.CreateLayer('line', geom_type=ogr.wkbLineString)
field_testfield = ogr.FieldDefn("fld_a", ogr.OFTString)
field_testfield.SetWidth(50)
layer.CreateField(field_testfield)
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField("fld_a",'test')
feature.SetGeometry(lineogr)
layer.CreateFeature(feature)
feature = None
ds = None

But it is easier with Fiona (another OGR Python bindings)

import fiona
import json
# shapefile schema
schema = {'geometry': 'LineString','properties': {'fld_a': 'str:50'}}
with fiona.open('line2.shp', 'w', 'ESRI Shapefile', schema) as layer:
      layer.write({'geometry': json.loads(line), 'properties': {'fld_a': 'test'}})

schema = {'geometry': 'Point','properties': {'fld_a': 'str:50'}}
with fiona.open('point.shp', 'w', 'ESRI Shapefile', schema) as layer:
      layer.write({'geometry': json.loads(point), 'properties': {'fld_a': 'test'}})

schema = {'geometry': 'Polygon','properties': {'fld_a': 'str:50'}}
with fiona.open('polygon.shp', 'w', 'ESRI Shapefile', schema) as layer:
      layer.write({'geometry': json.loads(poly), 'properties': {'fld_a': 'test'}})

Thanks @gene - I like the Fiona way and am working on it - but when you write json.loads(line) do you actually mean json.loads(lineogr)? — auslander, Jul 10 '17 at 15:33
No, Fiona uses the GeoJSON format and dictionaries. json.loads(line)transform the string to a Python dictionary — gene, Jul 10 '17 at 15:41
Thanks gene, I've edited the results section of my original post, see the update for new errors using the Fiona approach. I'm sure I'm making a simple logical error somewhere but can't discern it. — auslander, Jul 10 '17 at 16:02
Now I have the script successfully writing into a file, but it only writes one row and overwrites the file over and over. Do you suggest creating the shapefile container in a separate step and then opening it using 'a' flag or is there a more elegant solution? — auslander, Jul 10 '17 at 17:29
I does not understand your problem. No problem for me: one geometry = one shapefile — gene, Jul 10 '17 at 17:36
I got it. I'm iterating over a json response, not casting a single geometry. When I nested the for loop inside of the with fiona.open() call, everything worked. Thanks again gene! — auslander, Jul 11 '17 at 13:35

How to create line/polygon shapefiles from geojson using gdal in Python?

1 Answers1