src/HYDROTools/changeCoords.py

   1
   2 import numpy as np
   3 import MEDLoader as ml
   4 import medcoupling as mc
   5
   6 from pyproj import CRS
   7 from pyproj import Transformer
   8
   9 def changeCoords(fileIn, fileOut, epsgIn=2154, epsgOut=2154, offsetXin=0, offsetYin=0, offsetXout=0, offsetYout=0):
  10     """
  11     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.
  12
  13     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)
  14
  15     When there is only an origin offset change, use the same code for epsgIn and epsgOut.
  16     The operation is: Xout = Xin + offsetXin -offsetXout, Yout = Yin + offsetYin -offsetYout.
  17     If the epsgIn and epsgOut codes are equal, the transformation of system of coordinates (pyproj) is not called.
  18
  19     parameters:
  20     fileIn: full path of input MED file
  21     fileOut: full path of output MED file
  22     epsgIn: integer code of the EPSG for the input system of coordinates (for instance, Lambert II is 27572)
  23     epsgOut: integer code of the EPSG for the output system of coordinates (for instance, Lambert 93 is 2154)
  24     offsetXin: X origin in the input system of coordinates (default 0)
  25     offsetYin: Y origin in the input system of coordinates (default 0)
  26     offsetXout: X origin in the output system of coordinates (default 0)
  27     offsetYout: Y origin in the output system of coordinates (default 0)
  28     return:
  29     O if OK
  30     """
  31     # TODO: do we need to transform the exception in return codes ?
  32     if epsgIn != epsgOut:
  33         crs_in = CRS.from_epsg(epsgIn)
  34         crs_out = CRS.from_epsg(epsgOut)
  35         transformer = Transformer.from_crs(crs_in, crs_out)
  36
  37     meshMEDFileRead = ml.MEDFileMesh.New(fileIn)
  38     coords = meshMEDFileRead.getCoords()
  39     nb_comp = coords.getNumberOfComponents()
  40
  41     vx = coords[:,0]
  42     vy = coords[:,1]
  43     avx = vx.toNumPyArray() + offsetXin
  44     avy = vy.toNumPyArray() + offsetYin
  45
  46     if nb_comp == 3:
  47         vz = coords[:,2]
  48         avz = vz.toNumPyArray()
  49
  50     if epsgIn != epsgOut:
  51         navx, navy = transformer.transform(avx, avy)
  52     else:
  53         navx = avx
  54         navy = avy
  55     navx = navx - offsetXout
  56     navy = navy - offsetYout
  57
  58     if nb_comp == 3:
  59         navxy = np.stack((navx, navy, avz), axis=-1)
  60         ncoords = mc.DataArrayDouble(navxy)
  61         ncoords.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
  62     else:
  63         navxy = np.stack((navx, navy), axis=-1)
  64         ncoords = mc.DataArrayDouble(navxy)
  65         ncoords.setInfoOnComponents(["X [m]", "Y [m]"])
  66     meshMEDFileRead.setCoords(ncoords)
  67
  68     meshMEDFileRead.write(fileOut, 2)
  69     return 0