Skip to Content

ellipse_wang.py

  1. import fiona
  2. import numpy as np
  3. from shapely.geometry import shape
  4. # extraction des point x,y du fichier shapefile
  5. x, y = zip(*[(shape(point['geometry']).x, shape(point['geometry']).y) for point in fiona.open('points.shp'])
  6. # calcul de la matrice de covariance de Wang et al. (2015)
  7. x = np.array(x) 
  8. y = np.array(y) 
  9. meanx = np.mean(x)
  10. meany = np.mean(y)
  11. n = len(x)
  12. # calul de la matrice (basée sur n)
  13. xd = x-meanx
  14. yd = y-meany
  15. xyw= sum(xd*yd)
  16. x2w = sum(np.power(xd,2))
  17. y2w = sum(np.power(yd,2))
  18. cov = np.array([[x2w, xyw],[xyw,y2w]])/len(x)
  19. print cov
  20. [[ 11394.1712941    2160.35998689]
  21.  [  2160.35998689   2441.06326772]]
  22. # extraction des valeurs et  vecteurs propres
  23. eigenval,eigenvec = np.linalg.eig(cov)
  24. # rayons de l'ellipse (racine carrée des valeurs proprez
  25. sigmay,sigmax = np.sqrt(eigenval)
  26. print sigmax, sigmay
  27. [ 44.12521625 109.03302185  ]
  28. # angle de l'ellipse à partir des vecteurs propres 
  29. theta = 90 - np.rad2deg(np.arccos(eigenvec[0, 0]))
  30. print theta
  31. 77.1191513516
  32. # aire de l'ellipse
  33. aire = np.pi * sigmax*sigmay
  34. print aire
  35. 16194.143806210201
  36. # excentricité de l'ellipse
  37. exc = np.sqrt(1 - ((min(sigmax,sigmay)**2)/(max(sigmax,sigmay)**2)))
  38. print exc
  39. 0.91445132916786864
  40.  
  41. # solution présentée par Wang et al(2015)
  42. B = np.sqrt(np.power(x2w - y2w, 2) + 4 * np.power(xyw, 2))
  43. sigma2 = np.sqrt(((x2w + y2w) + B) / (2 * n))
  44. sigma1= np.sqrt(((x2w + y2w) - B) / (2 * n))
  45. print sigma1,sigma2
  46. 44.1252162467 109.033021847
  47. # angles de l'ellipse
  48. print np.degrees(np.arctan(((x2w - y2w) + B)/(2*xyw))), np.degrees(np.arctan(((x2w - y2w) - B)/(2*xyw)))
  49. 77.1191513516 -12.8808486484
  • 1238 lectures


Theme by Dr. Radut.