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.

asked Jul 6, 2017 at 14:31
9
  • why not use the gdal/ogr geojson driver to read the json in? Commented Jul 6, 2017 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." Commented Jul 6, 2017 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 Commented Jul 6, 2017 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 Commented Jul 6, 2017 at 15:32
  • try geom.ExportToWkt() Commented Jul 6, 2017 at 15:34

1 Answer 1

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'}}) 
answered Jul 7, 2017 at 14:29
6
  • 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)? Commented Jul 10, 2017 at 15:33
  • 1
    No, Fiona uses the GeoJSON format and dictionaries. json.loads(line)transform the string to a Python dictionary Commented Jul 10, 2017 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. Commented Jul 10, 2017 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? Commented Jul 10, 2017 at 17:29
  • I does not understand your problem. No problem for me: one geometry = one shapefile Commented Jul 10, 2017 at 17:36

Your Answer

Draft saved
Draft discarded

Sign up or log in

Sign up using Google
Sign up using Email and Password

Post as a guest

Required, but never shown

Post as a guest

Required, but never shown

By clicking "Post Your Answer", you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.