src/HYDROTools/changeCoords.py

   1
   2 import shutil
   3 import numpy as np
   4 import MEDLoader as ml
   5 import medcoupling as mc
   6
   7 import pyproj
   8
   9
  10 def transform(epsgIn, epsgOut, x, y):
  11     """
  12     Converting projection
  13     """
  14     # In older version of pyproj CRS is not available
  15     has_CRS = pyproj.__version__ >= "2.0.0"
  16     if has_CRS:
  17         crs_in = CRS.from_epsg(epsgIn)
  18         crs_out = CRS.from_epsg(epsgOut)
  19         new_x, new_y = Transformer.from_crs(crs_in, crs_out).transform(x, y)
  20     else:
  21         crs_in = pyproj.Proj(init="epsg:{}".format(epsgIn))
  22         crs_out = pyproj.Proj(init="epsg:{}".format(epsgOut))
  23         new_x, new_y = pyproj.transform(crs_in, crs_out, x, y)
  24
  25
  26 def changeCoords(fileIn, fileOut, epsgIn=2154, epsgOut=2154, offsetXin=0, offsetYin=0, offsetXout=0, offsetYout=0):
  27     """
  28     change the coordinates of a MED mesh file, with an optional offset reference in Input, and an optional change of system of coordinates (for instance Lambert II --> Lambert 93). The system of coordinates are given with EPSG (European Petroleum Survey Group) integer codes. The transformation uses pyproj.
  29
  30     When there is a change of system of coordinates, input offset should be given in order to get the correct input coordinates before transformation. example, with an input file in LambertII and an origin offset of (500000, 6000000), to obtain the coordinates in Lambert93 without origin offset: changeCoords("medfileLambertII.med", "medfileLambert93.med", 27572, 2154, 500000, 6000000, 0, 0)
  31
  32     When there is only an origin offset change, use the same code for epsgIn and epsgOut.
  33     The operation is: Xout = Xin + offsetXin -offsetXout, Yout = Yin + offsetYin -offsetYout.
  34     If the epsgIn and epsgOut codes are equal, the transformation of system of coordinates (pyproj) is not called.
  35
  36     parameters:
  37     fileIn: full path of input MED file
  38     fileOut: full path of output MED file
  39     epsgIn: integer code of the EPSG for the input system of coordinates (for instance, Lambert II is 27572)
  40     epsgOut: integer code of the EPSG for the output system of coordinates (for instance, Lambert 93 is 2154)
  41     offsetXin: X origin in the input system of coordinates (default 0)
  42     offsetYin: Y origin in the input system of coordinates (default 0)
  43     offsetXout: X origin in the output system of coordinates (default 0)
  44     offsetYout: Y origin in the output system of coordinates (default 0)
  45     return:
  46     O if OK
  47     """
  48     # TODO: do we need to transform the exception in return codes ?
  49     if fileOut != fileIn:
  50         shutil.copyfile(fileIn, fileOut)
  51
  52
  53     meshMEDFileRead = ml.MEDFileMesh.New(fileIn)
  54     coords = meshMEDFileRead.getCoords()
  55     nb_comp = coords.getNumberOfComponents()
  56
  57     vx = coords[:,0]
  58     vy = coords[:,1]
  59     avx = vx.toNumPyArray() + offsetXin
  60     avy = vy.toNumPyArray() + offsetYin
  61
  62     if nb_comp == 3:
  63         vz = coords[:,2]
  64         avz = vz.toNumPyArray()
  65
  66     if epsgIn != epsgOut:
  67         navx, navy = transform(epsgIn, epsgOut, avx, avy)
  68     else:
  69         navx = avx
  70         navy = avy
  71     navx = navx - offsetXout
  72     navy = navy - offsetYout
  73
  74     print(navx, navy)
  75
  76     if nb_comp == 3:
  77         navxy = np.stack((navx, navy, avz), axis=-1)
  78         ncoords = mc.DataArrayDouble(navxy)
  79         ncoords.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
  80     else:
  81         navxy = np.stack((navx, navy), axis=-1)
  82         ncoords = mc.DataArrayDouble(navxy)
  83         ncoords.setInfoOnComponents(["X [m]", "Y [m]"])
  84     meshMEDFileRead.setCoords(ncoords)
  85
  86     meshMEDFileRead.write(fileOut, 0)
  87     return 0