Changeset 2267 in lmdz_wrf

Timestamp:

Dec 26, 2018, 7:38:41 PM (6 years ago)

Author:

lfita

Message:

Getting there...

Location:

trunk/tools

Files:

• module_scientific.f90 (modified) (1 diff)
• nc_var_tools.py (modified) (10 diffs)

trunk/tools/module_scientific.f90

@@ -4932 +4932 @@
   END SUBROUTINE spacepercen_within_reg
 
+  SUBROUTINE spacepercen(dxA, dyA, xCAvals, yCAvals, NAvertexmax, xBAvals, yBAvals, dxB, dyB, xCBvals,&
+    yCBvals, NBvertexmax, xBBvals, yBBvals, dxyB, strict, Ngridsin, gridsin, percentages)
+  ! Subroutine to compute the space-percentages of a series of grid cells (B) into another series of 
+  !   grid-cells (A)
+  ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x 
+  !   and y axes. 
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: dxA, dyA, NAvertexAma, dxB, dyB,          &
+      NBvertexmax, dxyB
+    REAL(r_k), DIMENSION(dxA,dyA), INTENT(in)            :: xCAvals, yCAvals
+    REAL(r_k), DIMENSION(dxB,dyB), INTENT(in)            :: xCBvals, yCBvals
+    REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals
+    REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals
+    LOGICAL, INTENT(in)                                  :: strict
+    INTEGER, DIMENSION(dxA,dyA), INTENT(out)             :: Ngridsin
+    INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out)      :: gridsin
+    REAL(r_k), DIMENSION(dxA,dyA,Ngridsin), INTENT(out)  :: percentages
+
+! Local
+   INTEGER                                               :: iv, ix, iy
+   INTEGER                                               :: Nvertex, NvertexAgrid, Ncoin, Nsort
+   CHARACTER(len=20)                                     :: DS
+   REAL(r_k)                                             :: areapoly, areagpoly, totarea, totpercent
+   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                :: vertexgrid, poinsin
+
+!!!!!!! Variables
+! dxA, dyA: shape of the matrix with the grid points A
+! xCAvals, yCAvals: coordinates of the center of the grid cells A
+! NAvertexmax: Maximum number of vertices of the grid cells A
+! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex)
+! dxB, dyB: shape of the matrix with the grid points B
+! xCBvals, yCBvals: coordinates of the center of the grid cells B
+! NBvertexmax: Maximum number of vertices of the grid cells B
+! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex)
+! strict: give an error if the area of the polygon is not fully covered
+! Ngridsin: number of grids from grid B with some extension within the grid cell A
+! gridsin: indices of B grids within the grids of A
+! percentages: percentages of area of cells B of each of the grids within the grid cell A
+
+    fname = 'spacepercen'
+
+    DO ix = 1, dimxA
+      DO iy = 1, dimyA
+
+      ! Getting grid vertices
+      Nvertex = 0
+      DO iv=1, NAvertexmax
+        IF (xBAvals(ix,iy,iv) /= fillvalI) THEN
+          Nvertex = Nvertex + 1
+        END IF
+      END DO
+      IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid)
+      ALLOCATE(vertexgrid(Nvertex,2))
+      vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex)
+      vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex)
+  
+      CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xBCvals, yBCvals, NBvertexmax, &
+        xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:))
+      IF (ALLOCATE(poinsin)) DEALLOCATE(poinsin)
+      ALLOCATE(poinsin(Ngridsin(ix,iy),2))
+
+      DO iv=1, Ngridsin(ix,iy)
+        poinsin(i,1) = gridsin(ix,iy,iv,1)
+        poinsin(i,2) = gridsin(ix,iy,iv,2)
+      END DO
+
+      CALL spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, yBBvals,       &
+        Ngridsin(ix,iy), poinsin, strict, percentages(ix,iy,:))
+
+      END DO
+    END DO
+
 END MODULE module_scientific

trunk/tools/nc_var_tools.py

@@ -27200 +27200 @@
             [dimxn]: name of the dimension along the x-axis
             [varxn]: name of the variable with the values of dimxn
+            [varxbndsn]: name of the variable with the bounds of the values of varxn
+            [dimyn]: name of the dimension along the y-axis
+            [varyn]: name of the variable with the values of dimyn
+            [varybndsn]: name of the variable with the bounds of the values of varyn
             [slicedim]: series of values to define the slicing of the variable
               * [integer]: which value of the dimension
@@ -27240 +27244 @@
 
     iif = 0
+    Nexp = 8
     for fileinf in files:
         expectargs = '[netcdfname]:[dimvals]'
@@ -27254 +27259 @@
                 quit(-1)
 
-            if len(dimvalsref) != 6:
+            if len(dimvalsref) != Nexp:
                 print errormsg
                 print '  ' + fname + ": wrong number of dimension-values for " +     \
                   "reference !!"
-                print '    6 values are expected and you provided', len(dimvalsref)
-                print "    values must be '[dimxn];[varxn];[slicedim];[dimyn];" +    \
-                  "[varyn];[slicedim]'"
+                print '    ', Nexp,' values are expected and you provided',          \
+                  len(dimvalsref)
+                print "    values must be '[dimxn];[varxn];[varxbndsn];[slicedim];"+ \
+                  "[dimyn];[varyn];[varybndsn];[slicedim]'"
                 quit(-1)
 
@@ -27267 +27273 @@
             refdimxn = dimvalsref[0]
             refvarxn = dimvalsref[1]
-            slicedimx = dimvalsref[2]
-            refdimyn = dimvalsref[3]
-            refvaryn = dimvalsref[4]
-            slicedimy = dimvalsref[5]
+            refvarxbndsn = dimvalsref[2]
+            slicedimx = dimvalsref[3]
+            refdimyn = dimvalsref[4]
+            refvaryn = dimvalsref[5]
+            refvarybndsn = dimvalsref[6]
+            slicedimy = dimvalsref[7]
 
             if not gen.searchInlist(oncref.dimensions, refdimxn):
@@ -27290 +27298 @@
                 quit(-1)
 
+            if not oncref.variables.has_key(refvarxbndsn):
+                print errormsg
+                print '  ' + fname + ": reference file '" + netcdfrefn + "' does " + \
+                  "not have variable dimension for x '" + refvarxbndsn + "' !!"
+                ncvar = list(oncref.varibales.keys)
+                ncvar.sort()
+                print '    available ones:', ncvar
+                quit(-1)
+
             if not gen.searchInlist(oncref.dimensions, refdimyn):
                 print errormsg
@@ -27308 +27325 @@
                 quit(-1)
 
+            if not oncref.variables.has_key(refvarybndsn):
+                print errormsg
+                print '  ' + fname + ": reference file '" + netcdfrefn + "' does " + \
+                  "not have variable dimension for x '" + refvarybndsn + "' !!"
+                ncvar = list(oncref.varibales.keys)
+                ncvar.sort()
+                print '    available ones:', ncvar
+                quit(-1)
 
             slicedict = {refdimxn: slicedimx, refdimyn: slicedimy}
             ovarx = oncref.variables[refvarxn]
+            ovarxbnds = oncref.variables[refvarxbndsn]
             ovary = oncref.variables[refvaryn]
+            ovarybnds = oncref.variables[refvarybndsn]
 
             refvarx, dims = SliceVarDict(ovarx,slicedict)
+            refvarbndsx, dims = SliceVarDict(ovarxbnds,slicedict)
             refvary, dims = SliceVarDict(ovary,slicedict)
+            refvarbndsy, dims = SliceVarDict(ovarybnds,slicedict)
+
+            reflon, reflat =  gen.lonlat2D(refvarx, refvary)
+            refdx = reflon.shape[1]
+            refdy = reflon.shape[0]
+            refmaxvert = refvarxbnds.shape[2]
 
         else:
@@ -27326 +27360 @@
                 quit(-1)
 
-            if len(dimvalsget) != 6:
+            if len(dimvalsget) != Nexp:
                 print errormsg
                 print '  ' + fname + ": wrong number of dimension-values for get !!"
-                print '    6 values are expected and you provided', len(dimvalsget)
-                print "    values must be '[dimxn];[varxn];[slicedim];[dimyn];" +    \
-                  "[varyn];[slicedim]'"
+                print '    ', Nexp,' values are expected and you provided',          \
+                  len(dimvalsget)
+                print "    values must be '[dimxn];[varxn];[varxbndsn];[slicedim];"+ \
+                  "[dimyn];[varyn];[varybndsn];[slicedim]'"
                 quit(-1)
 
@@ -27338 +27373 @@
             getdimxn = dimvalsget[0]
             getvarxn = dimvalsget[1]
-            slicedimx = dimvalsget[2]
-            getdimyn = dimvalsget[3]
-            getvaryn = dimvalsget[4]
-            slicedimy = dimvalsget[5]
+            getvarxbndsn = dimvalsget[2]
+            slicedimx = dimvalsget[3]
+            getdimyn = dimvalsget[4]
+            getvaryn = dimvalsget[5]
+            getvarybndsn = dimvalsget[6]
+            slicedimy = dimvalsget[7]
 
             if not gen.searchInlist(oncget.dimensions, getdimxn):
@@ -27361 +27398 @@
                 quit(-1)
 
+            if not oncget.variables.has_key(getvarxbndsn):
+                print errormsg
+                print '  ' + fname + ": get file '" + netcdfgetn + "' does " +       \
+                  "not have variable dimension for x '" + getvarxbndsn + "' !!"
+                ncvar = list(oncget.varibales.keys)
+                ncvar.sort()
+                print '    available ones:', ncvar
+                quit(-1)
+
             if not gen.searchInlist(oncget.dimensions, getdimyn):
                 print errormsg
@@ -27379 +27425 @@
                 quit(-1)
 
+            if not oncget.variables.has_key(getvarybndsn):
+                print errormsg
+                print '  ' + fname + ": get file '" + netcdfgetn + "' does " +       \
+                  "not have variable dimension for x '" + getvarybndsn + "' !!"
+                ncvar = list(oncget.varibales.keys)
+                ncvar.sort()
+                print '    available ones:', ncvar
+                quit(-1)
 
             slicedict = {getdimxn: slicedimx, getdimyn: slicedimy}
             ovarx = oncget.variables[getvarxn]
+            ovarxbnds = oncget.variables[getvarxbndsn]
             ovary = oncget.variables[getvaryn]
 
             getvarx, dims = SliceVarDict(ovarx,slicedict)
+            getvarxbnds, dims = SliceVarDict(ovarxbnds,slicedict)
             getvary, dims = SliceVarDict(ovary,slicedict)
+            getvarybnds, dims = SliceVarDict(ovarybnds,slicedict)
+
+            getlon, getlat =  gen.lonlat2D(getvarx, getvary)
+            getdx = getlon.shape[1]
+            getdy = getlon.shape[0]
+            getmaxvert = getvarxbnds.shape[2]
 
         iif = iif + 1
+    # Getting the space weights
+    reflont = reflon.transpose()
+    reflatt = reflat.transpose()
+    refvarxbndst = refvarxbnds.transpose()
+    refvarybndst = refvarybnds.transpose()
+    getlont = getlon.transpose()
+    getlatt = getlat.transpose()
+    getvarxbndst = getvarxbnds.transpose()
+    getvarybndst = getvarybnds.transpose()
+
+    Ngridsint, gridsint, percenst = fsci.module_scientific.spacepercen(dxa=refdx,    \
+      dya=refdy, xcavals=reflont, ycavals=reflatt, navertmax=refmaxvert,             \
+      xbAvals=refvarxbndst, ybAvals=refvarybndst,                                    \
+      dxb=getdx, dyb=getdy, xcbvals=getlont, ycbvals=getlatt, nbvertmax=getmaxvert,  \
+      xbbvals=getvarxbndst, ybbvals=getvarybndst, strict=strict)
+
+    Ngridsin = Ngridint.transpose()
+    gridsin = gridint.transpose()
+    percens = percenst.transpose()
+
+    # Opening new file
+    ofilen = fname + '.nc'
+
+    onewnc = NetCDFFile(ofilen, 'w')
+    # dimensions
+    newdim = onewnc.createDimension('lon',refdx)
+    newdim = onewnc.createDimension('lat',refdy)
+    newdim = onewnc.createDimension('bnds',refmaxvert)
+    newdim = onewnc.createDimension('gridin',np.max(gridsin))
+    newdim = onewnc.createDimension('coords',2)
+
+    # variable-dimensions
+    newvar = onewnc.createVariable('lon','f8',('lat','lon'))
+    newvar[:] = reflon
+    basicvardef(newvar,'longitude','Longitude','degrees_east')
+    newavr.setncattr('bounds','lon_bnds')
+
+    newvar = onewnc.createVariable('lon_bnds','f8',('bnds', 'lat','lon'))
+    newvar[:] = refvarxbnds
+    basicvardef(newvar,'longitude_bnds','Bounds of Longitude','degrees_east')
+
+    newvar = onewnc.createVariable('lat','f8',('lat','lon'))
+    newvar[:] = reflat
+    basicvardef(newvar,'latitude','Latitude','degrees_north')
+    newavr.setncattr('bounds','lat_bnds')
+
+    newvar = onewnc.createVariable('lat_bnds','f8',('bnds', 'lat','lon'))
+    newvar[:] = refvarybnds
+    basicvardef(newvar,'latitude_bnds','Bounds of Latitude','degrees_north')
+
+    # variables space-weight
+    newvar = onewnc.createVariable('Ngrid','i',('gridin', 'lat','lon'))
+    newvar[:] = Ngridsin
+    basicvardef(newvar, 'Ngrid', "number of grids cells grom 'get' laying within" +  \
+      " 'ref'",'-')
+    newavr.setncattr('coordinates','lon lat')
+
+    newvar = onewnc.createVariable('gridin','i',('coord', 'gridin', 'lat', 'lon'))
+    newvar[:] = Ngridsin
+    basicvardef(newvar,'gridin',"coordinates of the grids cells grom 'get' laying "+ \
+      "within 'ref'",'-')
+    newavr.setncattr('coordinates','lon lat')
+
+    newvar = onewnc.createVariable('gridpercen','i',('coord', 'gridin', 'lat', 'lon'))
+    newvar[:] = percens
+    basicvardef(newvar,'gridpercen',"percentages of the grids cells grom 'get' " +   \
+      "laying within 'ref'", '1')
+    newavr.setncattr('coordinates','lon lat')
 
     # Getting values
 
+    add_global_PyNCplot(onewnc, 'nc_var_tools', fname, '1.0')
+
     return
 
-fvals= 'reference_data.nc:lon;lon;-1;lat;lat;-1,get_data.nc:lon;lon;-1;lat;lat;-1'
+fvals= 'reference_data.nc:lon;lon;lon_bnds;-1;lat;lat;lat_bnds;-1,get_data.nc:lon;lon;lon_bnds;-1;lat;lat;lat_bnds;-1'
 area_weighted('yes:mean,min,max',fvals,'ct_values')