Skip to Content

traitement avec osgeo.gdal

  1. from osgeo import gdal
  2. from numpy import *
  3. m = 12609
  4. n= 8810
  5. fp = [[1,m,m,1],[1,1,n,n]]
  6. tp = [[1616411.91,1617482.52,1617503.82,1616433.22],[9173672.33,9173641.84,9174389.85,9174420.34]]
  7. pix = zip(fp[0],fp[1])
  8. coor = zip(tp[0],tp[1])
  9. print pix
  10. ([(1, 1), (12609, 1), (12609, 8810), (1, 8810)]
  11. print coor
  12. [(1616411.9099999999, 9173672.3300000001), (1617482.52, 9173641.8399999999), (1617503.8200000001, 9174389.8499999996), (1616433.22, 9174420.3399999999)
  13. ])
  14.  
  15. gcps = []
  16. # calcule des paramètres avec gdal.GCP (points d'appui)
  17. for index, txt in enumerate(pix):
  18. gcps[index].GCPPixel = pix[index][0]
  19. gcps[index].GCPLine = 8810-int(pix[index][1])
  20. gcps[index].GCPX = coor[index][0]
  21. gcps[index].GCPY = coor[index][1]
  22. gcps.append(gdal.GCP())
  23.  
  24. geotransform = gdal.GCPsToGeoTransform( gcps )
  25. print geotransform
  26. (1616433.1325852636, 0.084914736675147789, -0.0024185492110491138, 9174420.3424183074, -0.0024183058378076717, -0.084914292200888988)
  27. >>> M = np.matrix([[geotransform[1],geotransform[4],geotransform[0]],[geotransform[2],geotransform[5],geotransform[3]]])
  28. print M
  29. matrix([[ 8.49147367e-02, -2.41830584e-03, 1.61643313e+06],
  30. [ -2.41854921e-03, -8.49142922e-02, 9.17442034e+06]])