Changeset 1173 in lmdz_wrf for trunk/tools


Ignore:
Timestamp:
Oct 11, 2016, 4:25:28 PM (9 years ago)
Author:
lfita
Message:

Adding function `reproject': Function to reproject values to another one
Using new subroutines in `module_ForInterpolate.F90':

  • `LlInterpolateProjection?': Subroutine which provides the indices for a given interpolation of a projection
  • `var2D_IntProj': Subroutine to interpolate a 2D variable
  • `var3D_IntProj': Subroutine to interpolate a 3D variable
  • `var4D_IntProj': Subroutine to interpolate a 4D variable
  • `var5D_IntProj': Subroutine to interpolate a 5D variable
Location:
trunk/tools
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/module_ForInterpolate.F90

    r900 r1173  
    11! Module to interpolate values from a giving projection and pressure interpolation
    22! To be included in a python
    3 ! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForInterpolate.F90  >& run_f2py.log
     3! Follow compilation instructions from Makefile
     4! Content
     5!  LlInterpolateProjection: Subroutine which provides the indices for a given interpolation of a projection
     6!  var2D_IntProj: Subroutine to interpolate a 2D variable
     7!  var3D_IntProj: Subroutine to interpolate a 3D variable
     8!  var4D_IntProj: Subroutine to interpolate a 4D variable
     9!  var5D_IntProj: Subroutine to interpolate a 5D variable
    410MODULE module_ForInterpolate
    511
    … …  
    714720END SUBROUTINE CoarseInterpolateExact
    715721
     722SUBROUTINE LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy,   &
     723  pdimx, pdimy)
     724! Subroutine which provides the indices for a given interpolation of a projection
     725
     726  USE module_generic
     727
     728  IMPLICIT NONE
     729
     730  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     731  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy
     732  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
     733  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
     734  CHARACTER(LEN=50), INTENT(in)                          :: intkind
     735  REAL(r_k), DIMENSION(3,16,pdimx,pdimy), INTENT(out)    :: outLlw
     736
     737! Local
     738  INTEGER                                                :: i,j,iv, ix, iy
     739  REAL(r_k)                                              :: mindiffLl, dist
     740  REAL(r_k), DIMENSION(idimx,idimy)                      :: difflonlat
     741  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
     742  INTEGER, DIMENSION(2)                                  :: iLl
     743  CHARACTER(LEN=50)                                      :: fname
     744
     745!!!!!!! Variables
     746! idimx, idimy: dimension length of the input projection
     747! pdimx, pdimy: dimension length of the target projection
     748! in[lon/lat]: longitudes and latitudes of the target projection
     749! proj[lon/lat]: longitudes and latitudes of the target projection
     750! intkind: kind of interpolation
     751!   'npp': nearest neighbourgh
     752!   'dis': weighted distance, grid-output for SW, NW, NE, SE
     753! outLlw: output interpolation result
     754!    for point pi,pj; up to 16 different values of
     755!       1st: i-index in input projection
     756!       2nd: j-index in input projection
     757!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
     758  fname = 'LlInterpolateProjection'
     759
     760  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
     761  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
     762
     763  ! Using case outside loop to be more efficient
     764  SELECT CASE(TRIM(intkind))
     765    CASE ('dis')
     766      DO i=1, pdimx
     767        DO j=1, pdimy     
     768          ! Find point
     769          difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.)
     770          mindiffLl = MINVAL(difflonlat)
     771          ! Getting the four surrounding values
     772          iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
     773          IF ( projlon(i,j) < inlonv(iLl(1),iLl(2)) ) THEN
     774            outLlw(1,1,i,j) = MAX(iLl(1)-1,1)
     775            outLlw(1,2,i,j) = MAX(iLl(1)-1,1)
     776            outLlw(1,3,i,j) = MIN(iLl(1),idimx)
     777            outLlw(1,4,i,j) = MIN(iLl(1),idimx)
     778          ELSE
     779            outLlw(1,1,i,j) = MIN(iLl(1),1)
     780            outLlw(1,2,i,j) = MIN(iLl(1),1)
     781            outLlw(1,3,i,j) = MAX(iLl(1)+1,idimx)
     782            outLlw(1,4,i,j) = MAX(iLl(1)+1,idimx)
     783          END IF
     784          IF ( projlat(i,j) < inlatv(iLl(1),iLl(2)) ) THEN
     785            outLlw(2,1,i,j) = MIN(iLl(1)-1,1)
     786            outLlw(2,2,i,j) = MAX(iLl(1),idimy)
     787            outLlw(2,3,i,j) = MAX(iLl(1),idimy)
     788            outLlw(2,4,i,j) = MIN(iLl(1)-1,1)
     789          ELSE
     790            outLlw(2,1,i,j) = MIN(iLl(1),1)
     791            outLlw(2,2,i,j) = MAX(iLl(1)+1,idimy)
     792            outLlw(2,3,i,j) = MAX(iLl(1)+1,idimy)
     793            outLlw(2,4,i,j) = MIN(iLl(1),1)
     794          END IF
     795          ! Computing distances
     796          DO iv=1,4
     797            ix = INT(outLlw(1,iv,i,j))
     798            iy = INT(outLlw(2,iv,i,j))
     799            dist = SQRT( (projlon(i,j)-inlonv(ix,iy))**2 + (projlat(i,j)-inlatv(ix,iy))**2 )
     800            IF ( dist /= 0.) THEN
     801              outLlw(3,iv,i,j) = 1./dist
     802            ELSE
     803              outLlw(3,iv,i,j) = 1.
     804            END IF
     805          END DO
     806          ! Normalizing
     807          outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j)))
     808        END DO
     809      END DO
     810    CASE ('npp')
     811      DO i=1, pdimx
     812        DO j=1, pdimy     
     813          ! Find point
     814          difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.)
     815          mindiffLl = MINVAL(difflonlat)
     816          outLlw(1:2,1,i,j) = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
     817          outLlw(3,1,i,j) = 1.
     818        END DO
     819      END DO
     820  END SELECT
     821
     822END SUBROUTINE LlInterpolateProjection
     823
     824SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
     825  pdimx, pdimy)
     826! Subroutine to interpolate a 2D variable
     827
     828  USE module_generic
     829
     830  IMPLICIT NONE
     831
     832  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     833  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy
     834  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
     835  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
     836  CHARACTER(LEN=50), INTENT(in)                          :: intkind
     837  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: var2Din
     838  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(out)         :: varout
     839
     840! Local
     841  INTEGER                                                :: i, j, iv, ix, iy
     842  REAL(r_k)                                              :: w
     843  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
     844  CHARACTER(LEN=50)                                      :: fname
     845
     846!!!!!!! Variables
     847! idimx, idimy: dimension length of the input projection
     848! pdimx, pdimy: dimension length of the target projection
     849! in[lon/lat]: longitudes and latitudes of the target projection
     850! proj[lon/lat]: longitudes and latitudes of the target projection
     851! intkind: kind of interpolation
     852!   'npp': nearest neighbourgh
     853!   'dis': weighted distance, grid-output for SW, NW, NE, SE
     854! outLlw: output interpolation result
     855!    for point pi,pj; up to 16 different values of
     856!       1st: i-index in input projection
     857!       2nd: j-index in input projection
     858!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
     859! var2Din: 2D variable to interpolate
     860! varout: variable interpolated on the target projection
     861  fname = 'var2D_IntProj'
     862
     863  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
     864    pdimy)
     865
     866  SELECT CASE (intkind)
     867    CASE('dis')
     868      DO i=1, pdimx
     869        DO j=1, pdimx
     870          varout(i,j) = 0.
     871          DO iv=1, 4
     872            ix = INT(outLlw(1,iv,i,j))
     873            iy = INT(outLlw(2,iv,i,j))
     874            w =  outLlw(3,iv,i,j)
     875            varout(i,j) = varout(i,j) + w*var2Din(ix,iy)
     876          END DO
     877        END DO
     878      END DO
     879    CASE('npp')
     880      DO i=1, pdimx
     881        DO j=1, pdimx
     882          ix = INT(outLlw(1,1,i,j))
     883          iy = INT(outLlw(2,1,i,j))
     884          varout(i,j) = var2Din(ix,iy)
     885        END DO
     886      END DO
     887  END SELECT
     888
     889END SUBROUTINE var2D_IntProj
     890
     891SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
     892  pdimx, pdimy, d3)
     893! Subroutine to interpolate a 3D variable
     894
     895  USE module_generic
     896
     897  IMPLICIT NONE
     898
     899  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     900  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy, d3
     901  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
     902  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
     903  CHARACTER(LEN=50), INTENT(in)                          :: intkind
     904  REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(in)       :: var3Din
     905  REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(out)      :: varout
     906
     907! Local
     908  INTEGER                                                :: i, j, k, iv, ix, iy
     909  REAL(r_k)                                              :: w
     910  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
     911  CHARACTER(LEN=50)                                      :: fname
     912
     913!!!!!!! Variables
     914! idimx, idimy: dimension length of the input projection
     915! pdimx, pdimy: dimension length of the target projection
     916! in[lon/lat]: longitudes and latitudes of the target projection
     917! proj[lon/lat]: longitudes and latitudes of the target projection
     918! intkind: kind of interpolation
     919!   'npp': nearest neighbourgh
     920!   'dis': weighted distance, grid-output for SW, NW, NE, SE
     921! outLlw: output interpolation result
     922!    for point pi,pj; up to 16 different values of
     923!       1st: i-index in input projection
     924!       2nd: j-index in input projection
     925!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
     926! var3Din: 3D variable to interpolate
     927! varout: variable interpolated on the target projection
     928  fname = 'var3D_IntProj'
     929
     930  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
     931    pdimy)
     932
     933  SELECT CASE (intkind)
     934    CASE('dis')
     935      DO i=1, pdimx
     936        DO j=1, pdimx
     937          DO k=1, d3
     938            varout(i,j,k) = 0.
     939            DO iv=1, 4
     940              ix = INT(outLlw(1,iv,i,j))
     941              iy = INT(outLlw(2,iv,i,j))
     942              w =  outLlw(3,iv,i,j)
     943              varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k)
     944            END DO
     945          END DO
     946        END DO
     947      END DO
     948    CASE('npp')
     949      DO i=1, pdimx
     950        DO j=1, pdimx
     951          DO k=1, d3
     952            ix = INT(outLlw(1,1,i,j))
     953            iy = INT(outLlw(2,1,i,j))
     954            varout(i,j,k) = var3Din(ix,iy,k)
     955          END DO
     956        END DO
     957      END DO
     958  END SELECT
     959
     960END SUBROUTINE var3D_IntProj
     961
     962SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
     963  pdimx, pdimy, d3, d4)
     964! Subroutine to interpolate a 4D variable
     965
     966  USE module_generic
     967
     968  IMPLICIT NONE
     969
     970  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     971  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy, d3, d4
     972  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
     973  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
     974  CHARACTER(LEN=50), INTENT(in)                          :: intkind
     975  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(in)    :: var4Din
     976  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(out)   :: varout
     977
     978! Local
     979  INTEGER                                                :: i, j, k, l, iv, ix, iy
     980  REAL(r_k)                                              :: w
     981  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
     982  CHARACTER(LEN=50)                                      :: fname
     983
     984!!!!!!! Variables
     985! idimx, idimy: dimension length of the input projection
     986! pdimx, pdimy: dimension length of the target projection
     987! in[lon/lat]: longitudes and latitudes of the target projection
     988! proj[lon/lat]: longitudes and latitudes of the target projection
     989! intkind: kind of interpolation
     990!   'npp': nearest neighbourgh
     991!   'dis': weighted distance, grid-output for SW, NW, NE, SE
     992! outLlw: output interpolation result
     993!    for point pi,pj; up to 16 different values of
     994!       1st: i-index in input projection
     995!       2nd: j-index in input projection
     996!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
     997! var4Din: 4D variable to interpolate
     998! varout: variable interpolated on the target projection
     999  fname = 'var4D_IntProj'
     1000
     1001  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
     1002    pdimy)
     1003
     1004  SELECT CASE (intkind)
     1005    CASE('dis')
     1006      DO i=1, pdimx
     1007        DO j=1, pdimx
     1008          DO k=1, d3
     1009            DO l=1, d4
     1010              varout(i,j,k,l) = 0.
     1011              DO iv=1, 4
     1012                ix = INT(outLlw(1,iv,i,j))
     1013                iy = INT(outLlw(2,iv,i,j))
     1014                w =  outLlw(3,iv,i,j)
     1015                varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l)
     1016              END DO
     1017            END DO
     1018          END DO
     1019        END DO
     1020      END DO
     1021    CASE('npp')
     1022      DO i=1, pdimx
     1023        DO j=1, pdimx
     1024          ix = INT(outLlw(1,1,i,j))
     1025          iy = INT(outLlw(2,1,i,j))
     1026          DO k=1, d3
     1027            DO l=1, d4
     1028              varout(i,j,k,l) = var4Din(ix,iy,k,l)
     1029            END DO
     1030          END DO
     1031        END DO
     1032      END DO
     1033  END SELECT
     1034
     1035END SUBROUTINE var4D_IntProj
     1036
     1037
     1038SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, varout, idimx, idimy,    &
     1039  pdimx, pdimy, d3, d4, d5)
     1040! Subroutine to interpolate a 5D variable
     1041
     1042  USE module_generic
     1043
     1044  IMPLICIT NONE
     1045
     1046  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     1047  INTEGER, INTENT(in)                                    :: idimx, idimy, pdimx, pdimy, d3, d4, d5
     1048  REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in)          :: projlon, projlat
     1049  REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
     1050  CHARACTER(LEN=50), INTENT(in)                          :: intkind
     1051  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(in) :: var5Din
     1052  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(out):: varout
     1053
     1054! Local
     1055  INTEGER                                                :: i, j, k, l, m, iv, ix, iy
     1056  REAL(r_k)                                              :: w
     1057  REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
     1058  CHARACTER(LEN=50)                                      :: fname
     1059
     1060!!!!!!! Variables
     1061! idimx, idimy: dimension length of the input projection
     1062! pdimx, pdimy: dimension length of the target projection
     1063! in[lon/lat]: longitudes and latitudes of the target projection
     1064! proj[lon/lat]: longitudes and latitudes of the target projection
     1065! intkind: kind of interpolation
     1066!   'npp': nearest neighbourgh
     1067!   'dis': weighted distance, grid-output for SW, NW, NE, SE
     1068! outLlw: output interpolation result
     1069!    for point pi,pj; up to 16 different values of
     1070!       1st: i-index in input projection
     1071!       2nd: j-index in input projection
     1072!       3rd: weight for i-index, j-index to use for ponderation (0<1.)
     1073! var5Din: 5D variable to interpolate
     1074! varout: variable interpolated on the target projection
     1075  fname = 'var5D_IntProj'
     1076
     1077  CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
     1078    pdimy)
     1079
     1080  SELECT CASE (intkind)
     1081    CASE('dis')
     1082      DO i=1, pdimx
     1083        DO j=1, pdimx
     1084          DO k=1, d3
     1085            DO l=1, d4
     1086              DO m=1, d5
     1087              varout(i,j,k,l,m) = 0.
     1088                DO iv=1, 4
     1089                  ix = INT(outLlw(1,iv,i,j))
     1090                  iy = INT(outLlw(2,iv,i,j))
     1091                  w =  outLlw(3,iv,i,j)
     1092                  varout(i,j,k,l,m) = varout(i,j,k,l,m) + w*var5Din(ix,iy,k,l,m)
     1093                END DO
     1094              END DO
     1095            END DO
     1096          END DO
     1097        END DO
     1098      END DO
     1099    CASE('npp')
     1100      DO i=1, pdimx
     1101        DO j=1, pdimx
     1102          ix = INT(outLlw(1,1,i,j))
     1103          iy = INT(outLlw(2,1,i,j))
     1104          DO k=1, d3
     1105            DO l=1, d4
     1106              DO m=1, d5
     1107                varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m)
     1108              END DO
     1109            END DO
     1110          END DO
     1111        END DO
     1112      END DO
     1113  END SELECT
     1114
     1115END SUBROUTINE var5D_IntProj
     1116
    7161117SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy,     &
    7171118  Ninpts)

  • trunk/tools/nc_var.py

    r1092 r1173  
    4444  'get_namelist_vars', 'get_Variables',                                              \
    4545  'getvalues_lonlat', 'grattr',                                                      \
    46   'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs', 'isvattrs', 'ivars',  \
    47   'ivattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics',                        \
     46  'grmattr', 'idims', 'igattrs', 'increaseDimvar', 'isgattrs',                       \
     47  'isvattrs', 'ivars', 'ivattrs', 'LMDZ_toCF', 'maskvar', 'model_characteristics',   \
    4848  'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation',    \
    49   'netcdf_fold_concatenation_HMT', 'Partialmap_Entiremap',                           \
     49  'netcdf_fold_concatenation_HMT', 'reproject', 'Partialmap_Entiremap',              \
    5050  'Partialmap_EntiremapFor', 'Partialmap_EntiremapForExact',                         \
    5151  'pinterp', 'remapnn',                                                              \
    … …  
    247247elif oper == 'Partialmap_EntiremapForExact':
    248248    ncvar.Partialmap_EntiremapForExact(opts.values, opts.ncfile, opts.varname)
     249elif oper == 'reproject':
     250    ncvar.reproject(opts.values, opts.ncfile, opts.varname)
    249251elif oper == 'seasmean':
    250252    ncvar.seasmean(timename, opts.ncfile, opts.varname)

  • trunk/tools/nc_var_tools.py

    r1170 r1173  
    9090# remapnn: Function to remap to the nearest neightbor a variable using projection from another file
    9191# remapnn_old: Function to remap to the nearest neightbor a variable using projection from another fil
     92# reproject: Function to reproject values to another one
    9293# seasmean: Function to compute the seasonal mean of a variable
    9394# sellonlatbox: Function to select a lotlan box from a data-set
    … …  
    1578215783#Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
    1578315784
     15785def reproject(values, filen, variable):
     15786    """ Function to reproject values to another one
     15787      reproject(values, filen, variable)
     15788      values=[inlonn],[inlatn],[projfile],[projlonn],[projlatn],[kind],[dimvars]
     15789        inlonn= name of the longitude values in the file to interpolate
     15790        inlatn= name of the latitude values in the file to interpolate
     15791        projlonn= name of the longitude values in the file with the target projection
     15792        projlatn= name of the latitude values in the file with the target projection
     15793        kind= kind of interpolation
     15794          'npp': nearest point
     15795          'dis': distance weighted (4 closest)
     15796        dimvars= ':', separated list of couples of dimension, and variable with its values [dimn]@[variabledim]
     15797      filen= name of file to interpolate
     15798      variable = ',' list of variable to interpolate ('all' for all variables)
     15799    """
     15800    import module_ForInt as fin
     15801    fname = 'reproject'
     15802
     15803    ofile = fname + '.nc'
     15804
     15805    if values == 'h':
     15806        print fname + '_____________________________________________________________'
     15807        print reproject.__doc__
     15808        quit()
     15809
     15810    arguments = '[inlonn],[inlatn],[projfile],[projlonn],[projlatn],[kind],[dimvars]'
     15811    gen.check_arguments(fname, values, arguments, ',')
     15812   
     15813    inlonn = values.split(',')[0]
     15814    inlatn = values.split(',')[1]
     15815    projfile = values.split(',')[2]
     15816    projlonn = values.split(',')[3]
     15817    projlatn = values.split(',')[4]
     15818    kind = values.split(',')[5]
     15819    dimvars = gen.Str_list(values.split(',')[6],':')
     15820
     15821    varns = gen.Str_list(variable,',')
     15822
     15823    inf = NetCDFFile(filen, 'r')
     15824    if not os.path.isfile(projfile):
     15825        print errormsg
     15826        print '  ' + fname + ": file with target projection '" + projfile +          \
     15827          "' does not exist !!"
     15828        quit(-1)
     15829    projf = NetCDFFile(projfile, 'r')
     15830
     15831    if not inf.variables.has_key(inlonn):
     15832        print errormsg
     15833        print '  ' +fname+ ": input file '" + filen + "' does not have longitude '" +\
     15834          inlonn + "' !!"
     15835        quit(-1)
     15836    if not inf.variables.has_key(inlatn):
     15837        print errormsg
     15838        print '  ' +fname+ ": input file '" + filen + "' does not have latitude '" + \
     15839          inlatn + "' !!"
     15840        quit(-1)
     15841    if not projf.variables.has_key(projlonn):
     15842        print errormsg
     15843        print '  ' + fname + ": file with target projection '" + projfile +          \
     15844          "' does not have longitude '" + projlonn + "' !!"
     15845        quit(-1)
     15846    if not projf.variables.has_key(projlatn):
     15847        print errormsg
     15848        print '  ' + fname + ": file with target projection '" + projfile +          \
     15849          "' does not have latitude '" + projlatn + "' !!"
     15850        quit(-1)
     15851
     15852    oinlon = inf.variables[inlonn]
     15853    oinlat = inf.variables[inlatn]
     15854    oprojlon = projf.variables[projlonn]
     15855    oprojlat = projf.variables[projlatn]
     15856
     15857    inlonvals0 = oinlon[:]
     15858    inlatvals0 = oinlat[:]
     15859    olonvals0 = oprojlon[:]
     15860    olatvals0 = oprojlat[:]
     15861
     15862    # Getting 2D longitudes and latitudes
     15863    ilonvals, ilatvals = gen.lonlat2D(ilonvals0, ilatvals0)
     15864    olonvals, olatvals = gen.lonlat2D(olonvals0, olatvals0)
     15865
     15866    if len(ilonvals.shape) != 2:
     15867        print errormsg
     15868        print '  ' + fname + ': wrong input longitude rank:',  ilonvals.shape,       \
     15869          'it must be 2!!'
     15870        quit(-1)
     15871
     15872    if len(ilatvals.shape) != 2:
     15873        print errormsg
     15874        print '  ' + fname + ': wrong input latitude rank:',  ilatvals.shape,        \
     15875          'it must be 2!!'
     15876        quit(-1)
     15877
     15878    if len(olonvals.shape) != 2:
     15879        print errormsg
     15880        print '  ' + fname + ': wrong output longitude rank:', olonvals.shape,       \
     15881          'it must be 2!!'
     15882        quit(-1)
     15883
     15884    if len(olatvals.shape) != 2:
     15885        print errormsg
     15886        print '  ' + fname + ': wrong output latitude rank:', olatvals.shape,        \
     15887          'it must be 2!!'
     15888        quit(-1)
     15889
     15890    for i in range(2):
     15891        if ilonvals[i] != ilatvals[i]:
     15892            print errormsg
     15893            print '  ' + fname + ': dimension',i,'of input longitude:', ilonvals[i], \
     15894              'differs from latitude:', ilatvals[i], 'they must be the same!!'
     15895            quit(-1)
     15896        if rlonvals[i] != rlatvals[i]:
     15897            print errormsg
     15898            print '  ' + fname + ': dimension',i,'of output longitude:', ilonvals[i],\
     15899              'differs from latitude:', ilatvals[i], 'they must be the same!!'
     15900            quit(-1)
     15901   
     15902    idx = ilonvals.shape[1]
     15903    idy = ilonvals.shape[0]
     15904    rdx = rlonvals.shape[1]
     15905    rdy = rlonvals.shape[0]
     15906
     15907    ilon = ilonvals.astype('float64')
     15908    ilat = ilatvals.astype('float64')
     15909    rlon = rlonvals.astype('float64')
     15910    rlat = rlatvals.astype('float64')
     15911
     15912    ilont = ilon.transpose()
     15913    ilatt = ilat.transpose()
     15914    rlont = rlon.transpose()
     15915    rlatt = rlat.transpose()
     15916
     15917    # dimension-vardims
     15918    vardimns = {}
     15919    for vardn in dimvars:
     15920        dimn = vardn.split('@')[0]
     15921        dimvn = vardn.split('@')[1]
     15922        vardimns[dimn] = dimvn
     15923
     15924    # File creation
     15925    onewnc = NetCDFFile(ofile, 'w')
     15926
     15927    # Dimensions
     15928    newdim = onewnc.createDimension('lon',rdx)
     15929    newdim = onewnc.createDimension('lat',rdy)
     15930
     15931    # Which variables to re-project
     15932    if varns[0] == 'all':
     15933        varns = inf.variables.keys()
     15934        # No interpolation of lon, lat
     15935        varns.remove(inlonn)
     15936        varns.remove(inlatn)
     15937
     15938    for varn in varns:
     15939        if not inf.variables.has_key(inlatn):
     15940            print errormsg
     15941            print '  ' + fname + ": input file '" + filen + "' does not have '" +    \
     15942              varn + "' !!"
     15943            quit(-1)
     15944        print '  ' + fname + ": re-projecting '" + varn + "' ..."
     15945        ovarin = inf.variables[varn]
     15946        varin = ovarin[:]
     15947
     15948        Nvarind = len(varin.shape)
     15949        if Nvarind < 2:
     15950            print errormsg
     15951            print '  ' + fname + ': wrong input data rank:', varin.shape,            \
     15952              'it must be at least 2!!'
     15953            quit(-1)
     15954
     15955        varin = varin.astype('float64')
     15956        varint = varin.transpose()
     15957
     15958        if Nvarind == 2:
     15959            varout = fin.var2D_IntProj(var2Din=varint, inlonv=inlont, inlatv=inlatt, \
     15960              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
     15961              pdimx=rdx, pdimy=rdy)
     15962        elif Nvarind == 3:
     15963            varout = fin.var3D_IntProj(var3Din=varint, inlonv=inlont, inlatv=inlatt, \
     15964              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
     15965              pdimx=rdx, pdimy=rdy, d3=varin.shape[0])
     15966        elif Nvarind == 4:
     15967            varout = fin.var4D_IntProj(var4Din=varint, inlonv=inlont, inlatv=inlatt, \
     15968              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
     15969              pdimx=rdx, pdimy=rdy, d3=varin.shape[0], d4=varin.shape[1])
     15970        elif Nvarind == 5:
     15971            varout = fin.var4D_IntProj(var4Din=varint, inlonv=inlont, inlatv=inlatt, \
     15972              projlon=rlont, projlat=rlatt, intkind=kind, idimx=idx, idimy=idy,      \
     15973              pdimx=rdx, pdimy=rdy, d3=varin.shape[0], d4=varin.shape[1],            \
     15974              d5=varin.shape[2])
     15975        else:
     15976            print errormsg
     15977            print '  ' + fname + ': rank:', Nvarind, 'for input data not ready !!'
     15978            quit(-1)
     15979   
     15980        # Variable's dimensions and their related variables
     15981        dimsfile = onewnc.dimensions()
     15982        for dimn in ovarin.dimensions():
     15983            if not gen.searchInlist(dimsfile, dimn):
     15984                odimn = inf.dimensions[dimn]
     15985                if odimn.isunlimited():
     15986                    newdim = onewnc.createDimension(dimn, None)
     15987                else:
     15988                    newdim = onewnc.createDimension(dimn, len(odimn))
     15989                if not vardimsn.has_key(dimn):
     15990                    print errormsg
     15991                    print '  ' + fname + ": dimension '" + dimn + "' does not have "+\
     15992                      'assigned any variable !!'
     15993                    print '    assigned pairs (via [dimvars]) _______'
     15994                    gen.print_dictionary(vardimsn)
     15995                    quit(-1)
     15996                else:
     15997                    add_vars(inf,onewnc,vardimsn[dimn])
     15998
     15999        # Interpolated variable
     16000        vardims = ovarin.dimensions()
     16001        newvardims = gen.replace_list(vardims,inlonn,'lon')
     16002        newvardims = gen.replace_list(newvardims,inlatn,'lat')
     16003        varattrs = ovarin.ncattrs()
     16004        if gen.searchInlist(varattrs,'_FillValue'):       
     16005            fillv = ovarin.getncattr('_FillValue')
     16006            newvar = onewnc.createVariable(varn, 'f4', tuple(newvardims),            \
     16007              fill_value=fillv)
     16008            varattrs.remove('_FillValue')
     16009        else:
     16010            newvar = onewnc.createVariable(varn, 'f4', tuple(newvardims))
     16011
     16012        newvar[:] = varout[:]
     16013        for attrn in varattrs:
     16014            attrv = ovarin.getncattr(attrn)
     16015            set_attribute(newvar, attrn, attrv)
     16016
     16017        onewnc.sync()
     16018
     16019    # Global attributes
     16020    add_globalattrs(inf,onewnc,'all')
     16021    onewnc.sync()
     16022
     16023    inf.close()
     16024    projf.close()
     16025    onewnc.close()
     16026
     16027    print fname + ": successful re-projection saved in '" + ofile + "' !!"
     16028
     16029    return
     16030
    1578416031def lonlatvarsFile(lonvals,latvals,dnx,dny,oLlnc):
    1578516032    """ Function to CF-define the longitudes and latitudes variables within a file
Note: See TracChangeset for help on using the changeset viewer.