Skip to Content

concave_hull1.py

  1. #!/usr/bin/env python
  2. # encoding: utf-8
  3. #title           :concave_hull1.py
  4. #description     :script accompagnant l'article  "Sur la création des enveloppes concaves (concave hull) et les divers moyens d'y parvenir".
  5. #auteur          :M.Laloux
  6. #date            :20140930 
  7. #==============================================================================
  8. from shapely.geometry import shape,mapping MultiLineString
  9. from shapely.ops import polygonize, unary_union
  10. from scipy.spatial import Delaunay
  11. import numpy as np
  12. import fiona
  13. couche = fiona.open("points_test.shp")
  14. points = [pt['geometry']['coordinates'] for pt in couche]
  15. points= np.array(points)
  16. # Delaunay
  17. tri = Delaunay(np.array(points))
  18.  
  19. # calcul des alpha_shapes
  20. def add_edge(edges, edge_points,i, j):
  21.      """Ajoute une ligne entre le i ème point et le j ème points si elle n'est pas dans la liste"""
  22.      if (i, j) in edges or (j, i) in edges:
  23.         # déjà dans la liste
  24.         return
  25.      edges.add( (i, j) )
  26.      edge_points.append(points[ [i, j] ])
  27.  
  28. def calc_alpha2(alpha):
  29.      edges = set()
  30.      edge_points = []
  31.      # loop over triangles:
  32.      for ia, ib, ic in tri.vertices:
  33.          # extraction des points de Delaunay
  34.          pa = points[ia]
  35.          pb = points[ib]
  36.          pc = points[ic]
  37.          # Longueurs des cotés du triangle
  38.          a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
  39.          b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
  40.          c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
  41.          # Semiperimètre du triangle
  42.          s = (a + b + c)/2.0
  43.          # Surface du triangle par la formule de Heron
  44.          area = np.sqrt(s*(s-a)*(s-b)*(s-c))
  45.          # rayon de filtrage
  46.          circum_r = a*b*c/(4.0*area)
  47.          if circum_r < alpha:
  48.             add_edge(edges,edge_points, ia, ib)
  49.             add_edge(edges,edge_points,ib, ic)
  50.             add_edge(edges,edge_points,ic, ia)
  51.      return edge_points
  52.  
  53. for alpha in np.arange(30,40):
  54.     m = MultiLineString(calc_alpha2(alpha))
  55.     triangles = list(polygonize(m))
  56.     a = unary_union(triangles)
  57.     # sélection de la première géométrie de type Polygon rencontrée
  58.     if a.geom_type == "Polygon":
  59.             print alpha
  60.             schema = {'geometry': 'Polygon','properties': {'test': 'str'}}
  61.             with fiona.open('myshp.shp', 'w', 'ESRI Shapefile', schema, crs) as layer:
  62.                 elem = {}
  63.                 elem['geometry'] = mapping(a) 
  64.                 elem['properties'] = {'test': '145'}
  65.                 layer.write(elem)
  66.             break
  67.     else:
  68.             continue
  • 1772 lectures


Theme by Dr. Radut.