• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# JGomezDans

a guest Jan 16th, 2009 471 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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top