2

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.

enter image description here

  • Do you want merge the shapefiles or union ? It is different – gene May 7 '17 at 12:25
0

Merging shapefiles and Unioning shapefiles is not the same thing.

1) Methods of Merging layers/shapefiles using OGR in python or Python GDAL/OGR Cookbook: Merge OGR Layers

2) Union multiple shapefiles with overlapping shapes

And it is easiest with Python modules as Shapely (unary_union) or GeoPandas (unary_union) (you cannot use ogr.UnionCascaded() because the command operates only on a multipolygon geometry)

With GeoPandas for example:

enter image description here

import geopandas as gpd
censuslayer = gpd.read_file("census.shp")
censuslayer
            geometry                                  id
0  LINESTRING (2.319786269430139 -233.45854922279...   1
1  LINESTRING (0.4904468911917945 -405.4164507772...   2
2  LINESTRING (1421.887143782383 -756.64961139896...   3
3  LINESTRING (2.319786269430139 -601.15576424870...   4
4  LINESTRING (626.1245142487047 -932.26619170984...   5
5  LINESTRING (844.3531060014478 -372.65400012735...   6

connectorLayer = gpd.read_file("connector.shp")
connectorlayer
            geometry                                  id
0  LINESTRING (996.3545156508733 -767.55209690988...  10
1  LINESTRING (108.2708235190839 -408.63377545266...  20
2  LINESTRING (51.71194948186536 -233.45854922279...  30
3  LINESTRING (77.2938756755466 -599.685684064271...  40
4  LINESTRING (870.5883607000272 -981.15071731242...  50
5  LINESTRING (1169.438309585492 -1111.5414507772...  60
6  LINESTRING (407.1494012015867 -680.53137397542...  70

Merge shapefiles: 13 features - LineString

 import pandas as pd
 merged =  pd.concat([ censuslayer,connectorLayer],ignore_index=True) 
 merged
            geometry                                  id
 0   LINESTRING (2.319786269430139 -233.45854922279...   1
 1   LINESTRING (0.4904468911917945 -405.4164507772...   2
 2   LINESTRING (1421.887143782383 -756.64961139896...   3
 3   LINESTRING (2.319786269430139 -601.15576424870...   4
 4   LINESTRING (626.1245142487047 -932.26619170984...   5
 5   LINESTRING (844.3531060014478 -372.65400012735...   6
 6   LINESTRING (996.3545156508733 -767.55209690988...  10
 7   LINESTRING (108.2708235190839 -408.63377545266...  20
 8   LINESTRING (51.71194948186536 -233.45854922279...  30
 9   LINESTRING (77.2938756755466 -599.685684064271...  40
 10  LINESTRING (870.5883607000272 -981.15071731242...  50
 11  LINESTRING (1169.438309585492 -1111.5414507772...  60
 12  LINESTRING (407.1494012015867 -680.53137397542...  70

Save resulting shapefile

merged.to_file('merged.shp')

Union resulting shapefile

union  = merged['geometry'].unary_union
union # the result is one feature 
<shapely.geometry.multilinestring.MultiLineString object at 0x10f507810>
  • Thank you - now i just need to figure out how to install geopandas. Am on MacOS and struggling as many have before me with installing pyproj, but this should do the the trick. To be clear, in my reading and understanding of union, that is what i wanted - most examples are of polygon unions - but the idea is to merge two layers inheriting the attributes of both and split geometries where needed (where the connectors intersect with the census streets, split the census street). Your response should do the trick though and thanks for the examples. – Cord May 7 '17 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.