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. OK, I Understand
 
Top