8

I have a .shp file, and I need to get the dual graph from this file (that is to say, each polygon centroid should become a vertex, and two vertices are linked if their associated polygons are adjacent).

I managed to get the centroids from the polygons of my file.

How can I add lines between the vertices that should indeed be connected?

| improve this question | |
  • Are you looking for an out of the box tool in QGIS or a script? – whyzar Apr 5 '18 at 13:41
  • At this point, every solution is a good solution to me. – Alexandre Apr 5 '18 at 14:10
  • I can totally understand that however, an out of the box tool in QGIS may not be available. If you want a script, you will have to at the very least show an attempt of creating a script that can be helped with. – whyzar Apr 5 '18 at 14:15
  • What form do you want the lines? As a set of line features with little link back to the original polygons? – Spacedman Apr 5 '18 at 14:52
  • Related to Can GRASS compute the line graph representation of an input graph? – underdark♦ Apr 5 '18 at 17:07
12

You can use Python only (versions 2.7.x and 3.x) without QGIS:

enter image description here

1) With Fiona, shapely and itertools

import fiona
from shapely.geometry import MuliPolygon, LineString, mapping
Multi = MultiPolygon([shape(poly['geometry']) for poly in fiona.open("polygons.shp")])
# creation of the dual graph shapefile
import itertools
# schema of the dual graph shapefile
schema = {'geometry': 'LineString','properties': {'test': 'int'}}
with fiona.open('dual_graph.shp','w','ESRI Shapefile', schema) as e:
   for poly1,poly2 in  itertools.combinations(Multi, 2):
       if poly1.touches(poly2):
           e.write({'geometry':mapping(LineString([poly1.centroid, poly2.centroid])), 'properties':{'test':1}})

enter image description here

2) with GeoPandas,shapely and itertools

import geopandas as gpd 
polygon = gpd.GeoDataFrame.from_file("polygons.shp")
geoms = polygon['geometry'].tolist()
dual = gpd.GeoDataFrame(gpd.GeoSeries([LineString([poly1.centroid, poly2.centroid]) for poly1,poly2 in  itertools.combinations(geoms, 2) if poly1.touches(poly2)]),columns=['geometry'])
dual.to_file("dual_graph2.shp")

3) I suppose you can do the same thing with PyQgis and itertools (Python 2.7 for 2.x versions and 3.x for 3.x versions)

| improve this answer | |
  • The first code works fine! I might just add that it should be "as e:" instead of "as e::", and that in my case (I didn't want edges between polygons sharing only one point) I added the condition " (poly1.intersection(poly2)).geom_type!="Point" " to the if loop. – Alexandre Apr 6 '18 at 8:31

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.