DM_crs_trafo.py
1 ## @package python.demo.DM_crs_trafo
2 #
3 # example script showing how to perform a crs transformation with pyDM
4 #
5 
6 from osgeo import osr
7 from opals import pyDM
8 
9 def getWkt(epsgCode):
10  crs = osr.SpatialReference()
11  crs.ImportFromEPSG(epsgCode)
12  return crs.ExportToWkt()
13 
14 
15 # create source wkt string epgs code
16 sourceCRS = getWkt(25833) # UTM zone 33N
17 targetCRS = getWkt(31256) # MGI Gauss-Krueger East
18 
19 # create trafo object
20 trafo = pyDM.Trafo(sourceCRS, targetCRS)
21 
22 
23 # simple coordinate transformation ---------------------
24 x = 529570.001
25 y = 5338590.000
26 z = 223.063
27 targetCoords = trafo.transform(x, y, z) # transform coordinates returned as tuple
28 # targetCoords = trafo.transform( (x, y, z) ) # alternatively coordinates can be provided as tuple
29 # targetCoords = trafo.transform( [x, y, z] ) # or list
30 
31 sourceCoords = trafo.transform(targetCoords[0], targetCoords[1], targetCoords[2], inv=True) # reverse transformation
32 
33 # check if trafo and reverse trafo lead to the original coordinates
34 eps = 0.001 # check to mm precision
35 assert( abs(sourceCoords[0]-x) < eps and abs(sourceCoords[1]-y) < eps and abs(sourceCoords[2]-z) < eps)
36 
37 
38 # geometry object transformation -----------------------
39 pt = pyDM.Point(x, y, z) # create point object
40 pyDM.transform(trafo, pt) # inplace transformation of point
41 assert( abs(targetCoords[0]-pt.x) < eps and abs(targetCoords[1]-pt.y) < eps and abs(targetCoords[2]-pt.z) < eps)
42 
43 lf = pyDM.PolylineFactory() # create polyline factory
44 lf.addPoint(x,y,z)
45 lf.addPoint(x+50,y,z)
46 line = lf.getPolyline() # get actual polyline object
47 pyDM.transform(trafo, line) # inplace transformation of polyline (multipart lines are also supported)
48 assert( abs(targetCoords[0]-line[0].x) < eps and abs(targetCoords[1]-line[0].y) < eps and abs(targetCoords[2]-line[0].z) < eps)
49 
50 pf = pyDM.PolygonFactory() # create polygon factory
51 pf.addPoint(x, y, z)
52 pf.addPoint(x+50, y, z)
53 pf.addPoint(x+50, y+100, z)
54 pf.addPoint(x, y+100, z)
55 pf.closePart()
56 polygon = pf.getPolygon()
57 pyDM.transform(trafo, polygon) # inplace transformation of polygon