The question is about Shapely and Fiona in pure Python without QGIS ("using command line and/or shapely/fiona").
A solution is
from shapely import shape, mapping
import fiona
# schema of the new shapefile
schema = {'geometry': 'Polygon','properties': {'area': 'float:13.3','id_populat': 'int','id_crime': 'int'}}
# creation of the new shapefile with the intersection
with fiona.open('intersection.shp', 'w',driver='ESRI Shapefile', schema=schema) as output:
for crim in fiona.open('crime_stat.shp'):
for popu in fiona.open('population.shp'):
if shape(crim['geometry']).intersects(shape(popu['geometry'])):
area = shape(crim['geometry']).intersection(shape(popu['geometry'])).area
prop = {'area': area, 'id_populat' : popu['id'],'id_crime': crim['id']}
output.write({'geometry':mapping(shape(crim['geometry']).intersection(shape(popu['geometry']))),'properties': prop})
The original two layers and the resulting layer
Part of the resulting layer table
You can use a spatial index (rtree here, look at GSE: Fastest way to join many points to many polygons in python and Using Rtree spatial indexing with OGR)
Another solution is to use GeoPandas (= Pandas + Fiona + Shapely)
import geopandas as gpd
g1 = gpd.GeoDataFrame.from_file("crime_stat.shp")
g2 = gpd.GeoDataFrame.from_file("population.shp")
data = []
for index, crim in g1.iterrows():
for index2, popu in g2.iterrows():
if crim['geometry'].intersects(popu['geometry']):
data.append({'geometry': crim['geometry'].intersection(popu['geometry']), 'crime_stat':crim['crime_stat'], 'Population': popu['Population'], 'area':crim['geometry'].intersection(popu['geometry']).area})
df = gpd.GeoDataFrame(data,columns=['geometry', 'crime_stat', 'Population','area'])
df.to_file('intersection.shp')
# control of the results in mi case, first values
df.head() # image from a Jupiter/IPython notebook
Update
You need to understand the definition of the spatial predicates. I use here the JTS Topology suite
As you can see there are only intersections and no crosses nor disjoint here. Some definitions from the Shapely manual
object.crosses(other): Returns True if the interior of the object intersects the interior of the other but does not contain it, and the dimension of the intersection is less than the dimension of the one or the other.
object.disjoint(other): Returns True if the boundary and interior of the object do not intersect at all with those of the other.
object.intersects(other): Returns True if the boundary and interior of the object intersect in any way with those of the other.
You can control it by a simple script (there are other solution but this one is the simplest)
i = 0
for index, crim in g1.iterrows():
for index2, popu in g2.iterrows():
if popu['geometry'].crosses(crim['geometry']):
i= i+1
print i
and the result is 0
Therefore, you only need intersects here.
Your script becomes
data = []
for index1, crim in g1.iterrows():
for index2, popu in g2.iterrows():
if popu['geometry'].intersects(crim['geometry']): # objects overlaps to partial extent, not contained
area_int = popu['geometry'].intersection(crim['geometry']).area
area_crim = crim['geometry'].area
area_popu = popu['geometry'].area #
# popu['properties'] is for Fiona, not for Pandas
popu_count = popu['PPL_CNT']
popu_frac = (area_int / area_popu) * popu_count#
# you must include the geometry, if not, it is a simple Pandas DataFrame and not a GeoDataframe
# Fiona does not accept a tuple as value of a field 'id': (index1, index2)
data.append({'geometry': crim['geometry'].intersection(popu['geometry']), 'id1': index1, 'id2':index2 ,'area_crim': area_crim,'area_pop': area_popu, 'area_inter': area_int, 'popu_frac': popu_frac} )
df = gpd.GeoDataFrame(data,columns=['geometry', 'id1','id2','area_crim', 'area_pop','area_inter'])
df.to_file('intersection.shp')
df.head()
Result:
EDIT: The following answer only works for the 1-D case. To extend it to 2-D, you will need to parameterize your links by the along-road distance, and replace the x-coordinate with the parameterized length. However, I'm fairly confident this is doable much simpler with Geopandas.
It would be too hard to give hints in the comments, so here's a script that should give you what you want. It's not written for efficiency--probably you could get geopandas to do what you want with some finegaling, but here ya go. It's also not written very generally, but that could be done if you have more than one attribute.
import geopandas as gpd
from shapely.geometry import LineString
import numpy as np
slgeom = [[(0,0),(7,0)],[(7,0),(13,0)],[(13,0),(15,0)],[(15,0),(19,0)]]
geoms = []
for s in slgeom:
geoms.append(LineString(s))
properti = ['a','b','c','d']
split_lines = gpd.GeoDataFrame(geometry=geoms)
split_lines['property'] = properti
olgeom = [[(0,0),(5,0)],[(5,0),(7,0), (10,0)],[(10,0),(13,0),(15,0)],[(15,0),(19,0)]]
geoms = []
for o in olgeom:
geoms.append(LineString(o))
susc = [1,2,3,4]
original_lines = gpd.GeoDataFrame(geometry=geoms)
original_lines['susc'] = susc
# Do split lines
xs1 = []
attrib1 = []
for g, a in zip(split_lines.geometry.values, split_lines.property.values):
x = g.coords.xy[0].tolist()
xs1.extend(x)
try:
attrib1[-1] = a
except:
pass
attrib1.extend([a for l in range(len(x))])
# Do originals
xs2 = []
attrib2 = []
for g, a in zip(original_lines.geometry.values, original_lines.susc.values):
x = g.coords.xy[0].tolist()
xs2.extend(x)
try:
attrib2[-1] = a
except:
pass
attrib2.extend([a for l in range(len(x))])
# Create numpy array for sorting
allxs = list(set(xs1 + xs2))
x_forsort = []
a1_forsort = []
a2_forsort = []
for x in allxs:
try:
idx = xs1.index(x)
a1_forsort.append(attrib1[idx])
except:
a1_forsort.append(None)
try:
idx = xs2.index(x)
a2_forsort.append(attrib2[idx])
except:
a2_forsort.append(None)
forsort = np.transpose(np.array([allxs, a1_forsort, a2_forsort]))
# Now sort based on x value (1st column)
sorteds = forsort[forsort[:,0].argsort()]
# Work through the sorted lists to create segments with the appropriate attributes
# Store results in a dictionary
output = dict()
output['geometry'] = []
output['attrib_1'] = []
output['attrib_2'] = []
for i in range(len(sorteds)-1):
# Store shapely linestring
output['geometry'].append(LineString([(sorteds[i,0],0),(sorteds[i+1,0],0)]))
# Store attributes
if i == 0:
output['attrib_1'].append(sorteds[i,1])
output['attrib_2'].append(sorteds[i,2])
else:
if sorteds[i,1] is None:
output['attrib_1'].append(output['attrib_1'][-1])
else:
output['attrib_1'].append(sorteds[i,1])
if sorteds[i,2] is None:
output['attrib_2'].append(output['attrib_2'][-1])
else:
output['attrib_2'].append(sorteds[i,2])
# Convert back to geopandas dataframe
out_gdf = gpd.GeoDataFrame(output)
out_gdf.crs = original_lines.crs
Result:
out_gdf
Out[185]:
attrib_1 attrib_2 geometry
0 a 1 LINESTRING (0 0, 5 0)
1 a 2 LINESTRING (5 0, 7 0)
2 b 2 LINESTRING (7 0, 10 0)
3 b 3 LINESTRING (10 0, 13 0)
4 c 3 LINESTRING (13 0, 15 0)
5 d 4 LINESTRING (15 0, 19 0)
Best Answer
1) The problem is that the intersection of two polygons is not always a polygon
The result is 2 Polygons and a LineString therefore the error
2) A solution is to convert all the geometries to the same type
Now, you can convert to a GeoDataFrame
3) Or to use GeoPandas Overlay
Export the result