domingo, 29 de enero de 2017

Intersecciones en la misma capa con GeoPandas

Cuando se desea efectuar la intersección entre rasgos de la misma capa es conveniente, primero, usar el método 'combinations' del módulo 'itertools' para excluir las repeticiones. La situación que se va a considerar se encuentra a continuación; donde la capa tiene 50 % de transparencia para facilitar la observación de los solapamientos.



El código siguiente produce las diferentes intersecciones, como GeoDataFrame, en una sola línea:


1
2
3
4
5
6
7
8
import geopandas as gpd
import itertools

pol12 = gpd.GeoDataFrame.from_file("pyqgis_data/polygon12.shp")

geoms = pol12['geometry'].tolist()

intersections = gpd.GeoDataFrame(gpd.GeoSeries(poly[0].intersection(poly[1]) for poly in  itertools.combinations(geoms, 2) if poly[0].intersects(poly[1])), columns=['geometry'])

Después de ejecutado en la Python Consola de QGIS se tiene que con 'matplotlib':


1
2
3
4
>>>from matplotlib import pyplot as plt
>>>plt.ion()
>>>intersections.plot()
<matplotlib.axes._subplots.AxesSubplot object at 0x854cda6c>

es posible corroborar que las intersecciones son precisamente las esperadas:


Con el método 'to_file' de GeoPandas es posible guardar el GeoDataFrame en cualquiera de los formatos vectoriales soportados:


1
2
3
>>>import fiona
>>>fiona.supported_drivers
{'FileGDB': 'raw', 'ESRI Shapefile': 'raw', 'OpenFileGDB': 'r', 'PCIDSK': 'r', 'AeronavFAA': 'r', 'SUA': 'r', 'GPSTrackMaker': 'raw', 'ARCGEN': 'r', 'PDS': 'r', 'DGN': 'raw', 'GeoJSON': 'rw', 'GPKG': 'rw', 'MapInfo File': 'raw', 'Idrisi': 'r', 'GPX': 'raw', 'DXF': 'raw', 'BNA': 'raw', 'SEGY': 'r', 'GMT': 'raw'}


No hay comentarios: