intersección de dos archivos de Python o de línea de comandos

Question

Soy nuevo trabajando con datos de geo y tratando de responder a una pregunta similar a la publicada aquí.

Utilizando la línea de comandos y/o torneadas/fiona, en la mayoría de manera eficiente y elegante.

He creado dos archivos con sus respectivos atributos numéricos:

1. las zonas con los datos de la población
2. los distritos de policía con el crimen estadísticas(número único)

Son independientes.

Me gustaría volver a calcular los tales que ambas estimaciones están presentes en cada uno de los shapefiles, por ejemplo. estimación de la delincuencia dentro del área pequeña (tomar los datos de 2 y mire todas las áreas en 1 que hay en ella, cortar esos que se cruzan y hacer de las matemáticas).

Yo tengo el de matemáticas descubierto y atrapado buscando elegantes formas y eficiente de los formatos de datos para hacer las intersecciones y cortar la se superpone a los límites apropiados.

Actualización tras la respuesta de @gen

La adaptación del código no encontrar ninguna de las intersecciones en los datos, el cual está estructurado de la siguiente manera (y sorprendente, ya que es a partir de fuentes oficiales):

1. polSHP_upd: cuarteles de policía a través de SA
2. sal_subSHP: subconjunto de pequeñas regiones con datos de la población en Ciudad del Cabo

Los archivos en vivo aquí:

https://github.com/AnnaMag/GIS-analyses

Lo que me estoy perdiendo?

polSHP_upd = 'crime/GIS-analyses/polPres_updated4crimes/polPres_updated4crimes.shp'
sal_subSHP = 'crime/GIS-analyses/sal_sub/sal_sub.shp'

import geopandas as gpd
g1 = gpd.GeoDataFrame.from_file(polSHP_upd)
g2 = gpd.GeoDataFrame.from_file(sal_subSHP)
data = []
# A: region of the police data with criminal record
# C: small area with population data
# we look for all small areas intersecting a given C_i, calculate the  fraction of inclusion, scale the
# population accordingly: area(A_j, where A_j crosses C_i)/area(A_j)*  popul(A_j)
for index, crim in g1.iterrows():
   for index2, popu in g2.iterrows():      
      if popu['geometry'].crosses(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 # there is a Shape_Area field in properties: to check precision
         popu_count = popu['properties']['PPL_CNT']
         popu_frac = (area_int / area_popu) * popu_count# fraction of the pop area contained inside the crim
         data.append({'id': (index1, index2) ,'area_crim': area_crim,'area_pop': area_popu, 'area_inter': area_int, 'popu_frac': popu_frac} )
         # data.append({'geometry': crim['geometry'].intersection(popu['geometry']), 'crime_stat':crim['crime_stat'], 'Population': popu['Population'], 'area':crim['geometry'].intersection(popu['geometry']).area})

      elif popu['geometry'].intersects(crim['geometry']): 
        pass
        #print("intersect")

      elif popu['geometry'].disjoint(crim['geometry']): 
        pass`

Answer 1

La pregunta es bien formada y Fiona en Python puro sin QGIS ("el uso de la línea de comandos y/o torneadas/fiona").

Una solución es

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})

El original de dos capas y la capa resultante

Parte de la capa resultante de la tabla

Usted puede utilizar un índice espacial (rtree aquí, miren GSE: la forma más Rápida de juntar muchos puntos a muchos polígonos en python y el Uso de Rtree espacial de la indización con OGR)

Otra solución es utilizar GeoPandas (= Pandas + Fiona + bien formada)

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

Actualización

Usted necesita entender la definición de la ordenación de los predicados. Utilizo aquí la STC Topología suite

Como usted puede ver no sólo son las intersecciones y no cruza ni discontinuo aquí. Algunas definiciones de la bien formada manual

objeto.cruces(otros): Devuelve True si el interior del objeto cruza el interior de la otra, pero no la contienen, y la dimensión de la intersección es menor que la dimensión de la una o la otra.
objeto.discontinuo(otros): Devuelve True si la frontera y en el interior del objeto no se cruzan con los de los otros.
objeto.cruza(otros): Devuelve True si la frontera y en el interior del objeto se cruzan en cualquier forma con los de los otros.

Se puede controlar con un simple script (hay otra solución, pero éste es el más simple)

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

y el resultado es 0

Por lo tanto, usted necesita solamente se cruza aquí.

Su guión se convierte en

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()

Resultado:

Answer 2

Usted puede hacer que en QGIS, sin 'bien formada' y 'fiona', mediante el uso de PyQGIS. Para un arreglo similar de shapefiles (ver la siguiente imagen) de la respuesta en el enlace:

Cómo calcular el tamaño de un área en particular por debajo de un búfer en QGIS

Este código:

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)

funciona de manera adecuada para la producción de una capa de memoria con la intersección de las funciones. Los atributos de la tabla incluye las áreas requeridas de cada polígono; como se puede observar en la siguiente imagen: