Guest User

JGomezDans

a guest
Jan 16th, 2009
518
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy
  2.  
  3. def ConvertCoords ( x, y, z):
  4.     from osgeo import osr
  5.     # Set up spatial reference systems
  6.     # using EPSG 28992
  7.     proj = osr.SpatialReference()
  8.     proj.ImportFromEPSG(28992)
  9.     proj.SetTOWGS84(565.237, 50.0087, 465.658, -0.406857, 0.350733, -1.87035, 4.0812 )
  10.     # lat/lon WGS84
  11.     latlong = osr.SpatialReference()
  12.     latlong.ImportFromProj4('+proj=latlong +datum=WGS84')
  13.     #Define transform, from EPSG28992 to LatLon/WGS84
  14.     transform = osr.CoordinateTransformation( proj, latlong  )
  15.     (lon, lat,z) = transform.TransformPoint( x, y, z)
  16.     return (lon,lat, z)
  17.  
  18. if __name__=="__main__":
  19.     from osgeo import osr
  20.     # Set up spatial reference systems
  21.     # using EPSG 28992
  22.     proj = osr.SpatialReference()
  23.     proj.ImportFromEPSG(28992)
  24.     proj.SetTOWGS84(565.237, 50.0087, 465.658, -0.406857, 0.350733, -1.87035, 4.0812 )
  25.     # lat/lon WGS84
  26.     latlong = osr.SpatialReference()
  27.     latlong.ImportFromProj4('+proj=latlong +datum=WGS84')
  28.  
  29.     from osgeo import ogr
  30.     polygon =numpy.array ([[132817.006604435708141, 550302.852720651309937, 0.],\
  31.               [131182.28895997320069, 558340.214472591876984, 0.],\
  32.               [132578.61028128489852, 558748.893883707583882, 0.],\
  33.               [136631.347774848196423, 553436.061539204441942, 0.],\
  34.               [136631.347774848196423, 553436.061539204441942, 0.],\
  35.               [132817.006604435708141, 550302.852720651309937, 0.]])
  36. @@  # We can use two methods: convert individual points, or create an OGR feature and convert that 
  37. @@  # Method 1: Use convert individual points
  38.     result = [ConvertCoords ( polygon[i, 0], polygon[i, 1], polygon[i,2] ) for i in xrange(polygon.shape[0])]
  39.     print result
  40.     # Method 2: Convert whole polygon...
  41.     wkt = "POLYGON(("
  42.     edges = ["%f %f %f,"%( polygon[i,0], polygon[i,1], polygon[i,2]) for i in xrange(polygon.shape[0])]
  43.     wkt = wkt+"".join (edges[:-1])+edges[-1].replace(",", "")+"))"
  44.     geom = ogr.CreateGeometryFromWkt ( wkt )
  45.     geom.AssignSpatialReference ( proj )
  46.     geom.TransformTo ( latlong)
  47.     print geom.ExportToWkt()   
RAW Paste Data