# JGomezDans

a guest
Jan 16th, 2009
869
0
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()