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.

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(censusstreetlayersrc.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.