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.

Answer 1

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'}})