2

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 6 '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 6 '17 at 14:59
  • See pcjericks.github.io/py-gdalogr-cookbook/… but use GeoJSON instead of ESRI Shapefile in the driver = ogr.GetDriverByName("ESRI Shapefile") line – Ian Turton Jul 6 '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 6 '17 at 15:32
  • try geom.ExportToWkt() – Ian Turton Jul 6 '17 at 15:34
4

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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