[GIS] How to create line/polygon shapefiles from geojson using gdal in Python

gdalogropen-source-gispython

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.

Best Answer

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

Related Solutions

[GIS] Why do Shapefiles and GeoJSON behave differently in GDAL Python

Gross, but:

del lyrBuffer
del ds  # flushes Buffer.geojson to disk
ds = GeoJSONDriver.Open('Buffer.geojson', 0)
lyrBuffer = ds.GetLayerByIndex(0)
print lyrBuffer.GetFeatureCount()

From what I can tell, OGRGeoJSONDataSource.CreateLayer() returns a OGRGeoJSONWriteLayer, not a OGRGeoJSONLayer (which is what you get when opening an existing file). The former doesn't implement GetFeatureCount() (or much other useful stuff), and the datasource doesn't (a) actually flush the file to disk until it's deleted (ignoring SyncToDisk())or (b) expose any of the private methods that do flush it or enable proper re-reading of the layer.

By contrast, the Shapefile driver uses a single layer class for both reading and writing, and changes are propagated much better when you're editing.

Normally Fiona is a nicer python-based interface to OGR, but in this case it still requires reloading the file:

import fiona

with fiona.open('test.shp', 'r') as source:
    with fiona.open('Buffer.geojson', 'w',
                    crs=source.crs,
                    driver="GeoJSON",
                    schema=source.schema) as sink:

        for f in source:
            sink.write(f)

with fiona.open('Buffer.geojson', 'r') as sink:
    print len(sink)

[GIS] create polygon adjacency matrix using python

Just figured out a way to do this in R, inspired by the link I posted in the comments (which uses outdated functionalities but the right packages); as an example, I'll use the "Southwesternmost" counties in Michigan (shapefile here--Allegan, Berrien, Cass, Kalamazoo, St. Joseph, and Van Buren counties)

visual of sample counties

So, ordering counties alphabetically, the (length>0, i.e., ignoring corners) adjacency matrix we expect is:

  A B C K S V
A 0 0 0 1 0 1
B 0 0 1 0 0 1
C 0 1 0 0 1 1
K 1 0 0 0 1 1
S 0 0 1 1 0 0
V 1 1 1 1 0 0

To get this in R, we can do the following:

library(maptools)
library(spdep)

counties<-readShapeSpatial("~/Desktop/test.shp")
adj_mat<-nb2mat(poly2nb(counties,queen=F,row.names=counties$NAME),style="B")

> adj_mat
           [,1] [,2] [,3] [,4] [,5] [,6]
Allegan       0    0    0    1    0    1
Berrien       0    0    1    0    0    1
Cass          0    1    0    0    1    1
Kalamazoo     1    0    0    0    1    1
St. Joseph    0    0    1    1    0    0
Van Buren     1    1    1    1    0    0

For some bells and whistles:

For some purposes, being a neighbor to yourself is meaningful:

diag(adj_mat)<-rep(1,nrow(adj_mat))

May also want to define column names in addition to row names:

colnames(adj_mat)<-rownames(adj_mat)

Best Answer

Related Solutions

[GIS] Why do Shapefiles and GeoJSON behave differently in GDAL Python

[GIS] create polygon adjacency matrix using python

Related Question