I am using the Python GDAL/OGR bindings. I have two linestring shapefile layers. They are the census TIGER streets and a new set of line segments that represent lines connecting census block centroids to those streets. I want to merge the two into one shapefile (so I can do some route analysis) where the connector line segments break the street map segments where they intersect. I can do this with QGIS but need to automate this over thousands of datasets.
EDIT:
I believe I have GDAL 2.1 on a MacOS (Sierra)
$ geos-config --version
3.5.0
$ gdal-config --version
2.1.0
$ gdal-config --dep-libs
-L/Users/me/anaconda3/lib -lgeos_c -lsqlite3 -lexpat -L/Users/me/anaconda3/lib -lxerces-c -lpthread -L/Users/me/anaconda3/lib -lhdf5 -L/Users/me/anaconda3/lib -lmfhdf -ldf -L/Users/me/anaconda3/lib -lgeotiff -L/Users/me/anaconda3/lib -ltiff -L/usr/lib -lpq -lz -lpthread -ldl -L/Users/me/anaconda3/lib -lcurl -lssl -lcrypto -lssl -lcrypto -lz -liconv -lxml2 -lz -lpthread -licucore -lm -L/Users/me/anaconda3/lib -lkea -L/Users/me/anaconda3/lib -lhdf5 -lhdf5_hl -lhdf5_cpp
The below results in no errors but with an empty Shapefile with 0 features.
Edit 2: I do get an error. Segmentation fault: 11 on the union operation which is not an uncommon problem. I have verified the shapefile created using my CreateShapeFromCSV in QGIS and all seems fine.
def CreateNewShapeLayer(shapeName):
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource(shapeName)
# create the spatial reference, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
# create the layer
layer = data_source.CreateLayer(shapeName.rsplit("/", 1)[-1], srs, ogr.wkbLineString)
print ("Created layer {}".format(shapeName))
return layer
def UnionBlockCentroidStreetLines():
censusstreetlayersrc = "/Users/cthomas/Development/Data/spatial/Network/streets/tl_2016_06000_roads_la_clipped.shp"
censusStreets = ogr.Open(censusstreetlayersrc, 0)
censuslayer = censusStreets.GetLayer()
connectorlayersrc = "/Users/cthomas/Development/Data/spatial/Network/streets/street_segment_block_centroid_connectors.csv"
# Convert my CSV connector layer to a Shapefile
CreateShapeFromCSV(connectorlayersrc)
connectorStreets = ogr.Open(connectorlayersrc.replace(".csv", ".shp"), 0)
connectorLayer = connectorStreets.GetLayer()
# Create a new empty shapefile to receive the results of my union
outputLayer = CreateNewShapeLayer(censusstreetlayersrc.rsplit("/", 1)[0] + "/la_streets_with_block_centroid_connectors.shp")
print("About to union the layers")
censuslayer.Union(connectorLayer, outputLayer)
outputLayer = None
censuslayer = None
connectorLayer = None
UnionBlockCentroidStreetLines()
I see I can do something feature by feature, but was hoping I could do it for the whole layer.
The two layers are represented in the image below.
