Intersecting two shapefiles from Python or command line?

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:

You can do that in QGIS, without 'shapely' and 'fiona', by using PyQGIS. For a similar arrangement of shapefiles (see next image) from the answer in your link:

How to calculate the size of a particular area below a buffer in QGIS

This code:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feats0 = [feat for feat in layers[0].getFeatures()]
feats1 = [feat for feat in layers[1].getFeatures()]

geom_intersec = [ feats0[0].geometry().intersection(feat.geometry()).exportToWkt()
                  for feat in feats1 ] 

geom_int_areas = [ feats0[0].geometry().intersection(feat.geometry()).area()
                   for feat in feats1 ] 

crs = layers[0].crs()

epsg = crs.postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&field=area&index=yes"

intersections = QgsVectorLayer(uri, 
                      'intersections', 
                      'memory')

QgsMapLayerRegistry.instance().addMapLayer(intersections)

prov = intersections.dataProvider()

n = len(geom_intersec)

feats = [ QgsFeature() for i in range(n) ]

for i, feat in enumerate(feats): 
    feat.setGeometry(QgsGeometry.fromWkt(geom_intersec[i]))
    feat.setAttributes([i, geom_int_areas[i]])

prov.addFeatures(feats)

it works adequately for producing a memory layer with the intersection features. The attributes table includes the required areas of each polygon; as it can be observed at next image: