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.
1 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'}})
-
Thanks @gene - I like the Fiona way and am working on it - but when you write
json.loads(line)
do you actually meanjson.loads(lineogr)
?auslander– auslander2017年07月10日 15:33:03 +00:00Commented Jul 10, 2017 at 15:33 -
1No, Fiona uses the GeoJSON format and dictionaries.
json.loads(line)
transform the string to a Python dictionarygene– gene2017年07月10日 15:41:48 +00:00Commented 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.auslander– auslander2017年07月10日 16:02:37 +00:00Commented 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?auslander– auslander2017年07月10日 17:29:54 +00:00Commented Jul 10, 2017 at 17:29 -
I does not understand your problem. No problem for me: one geometry = one shapefilegene– gene2017年07月10日 17:36:55 +00:00Commented Jul 10, 2017 at 17:36
dataSource = driver.Open(jsonObj, 0)
and then didfor feature in layer: geom = feature.GetGeometryRef() \ print geom
it yieldedRuntimeError: not a string