Changeset 1620 in lmdz_wrf

Timestamp:

Sep 7, 2017, 2:52:14 PM (8 years ago)

Author:

lfita

Message:

Adding in `polygon_properties' calculation of:

• area-weighted center
• coordinates of the minimum and maximum values
• coordinates of the area-weighted and distance to minium/maximum center of the polygon

File:

• trunk/tools/module_scientific.f90 (modified) (13 diffs)

trunk/tools/module_scientific.f90

@@ -32 +32 @@
   CONTAINS
 
-  SUBROUTINE all_polygons_properties(dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN,       &
-    Npolyptss, xxtrms, yxtrms, meanctrs, areas, nvals, xvals, mvals, m2vals, stdvals, Nquant, quants)
+  SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN,  &
+    Npolyptss, xxtrms, yxtrms, meanctrs, meanwctrs, areas, nvals, xvals, mvals, m2vals, stdvals,      &
+    Nquant, quants, nvcoords, xvcoords, meanvnctrs, meanvxctrs)
 ! Subroutine to determine the properties of all polygons in a 2D field:
 !   Number of grid points
 !   grid-point coordinates of the minimum and maximum of the path along x,y axes
-!   center from the mean of the coordinates of the paths locations
+!   grid coordinates of center from the mean of the coordinates of the poygon locations
+!   lon, lat center from the area weighted mean of the coordinates of the polygon locations
 !   area of the polygon (km2)
 !   minimum and maximum of the values within the polygon
@@ -45 +47 @@
 !   number of quantiles
 !   quantiles of the values within the polygon
+!   grid coordinates of the minimum, maximum value within the polygon
+!   lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
 
   IMPLICIT NONE
 
+  LOGICAL, INTENT(in)                                    :: dbg
   INTEGER, INTENT(in)                                    :: dx, dy, Npoly, Nquant
   INTEGER, DIMENSION(dx,dy), INTENT(in)                  :: polys
   REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: lon, lat, values
   REAL(r_k), INTENT(in)                                  :: xres, yres
-  CHARACTER(len=50), INTENT(in)                          :: projN
+  CHARACTER(len=1000), INTENT(in)                        :: projN
   INTEGER, DIMENSION(Npoly), INTENT(out)                 :: Npolyptss
-  INTEGER, DIMENSION(Npoly, 2), INTENT(out)              :: xxtrms, yxtrms, meanctrs
+  INTEGER, DIMENSION(Npoly,2), INTENT(out)               :: xxtrms, yxtrms, meanctrs
   REAL(r_k), DIMENSION(Npoly), INTENT(out)               :: areas
   REAL(r_k), DIMENSION(Npoly), INTENT(out)               :: nvals, xvals, mvals, m2vals, stdvals
   REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out)       :: quants
+  INTEGER, DIMENSION(Npoly,2), INTENT(out)               :: nvcoords, xvcoords
+  REAL(r_k), DIMENSION(Npoly,2), INTENT(out)             :: meanwctrs, meanvnctrs, meanvxctrs
 
 ! Local
@@ -72 +79 @@
 ! projN: name of the projection
 !   'lon/lat': for regular longitude-latitude
+!   'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
 ! Npolyptss: number of points of the polygons
 ! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons
 ! meanctrs: center from the mean of the coordinates of the polygons
+! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons
 ! areas: area of the polygons [km]
 ! [n/x]vals: minimum and maximum of the values within the polygons
@@ -82 +91 @@
 ! Nquant: number of quantiles
 ! quants: quantiles of the values within the polygons
+! [n/x]vcoords: grid coordinates of the minimum/maximum of the values within the polygons
+! meanv[n/x]ctrs: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygons
 
   fname = 'all_polygons_properties'
@@ -90 +101 @@
   yxtrms = fillval64
   meanctrs = fillval64
+  meanwctrs = fillval64
   areas = fillval64
   nvals = -1
@@ -97 +109 @@
   stdvals = fillval64
   quants = fillval64
+  nvcoords = -1
+  xvcoords = -1
+  meanvnctrs = fillval64
+  meanvxctrs = fillval64
 
   DO ip=1, Npoly
     boolpoly = polys == ip
-    CALL polygon_properties(dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),     &
-      xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), areas(ip), nvals(ip), xvals(ip), mvals(ip),         &
-      m2vals(ip), stdvals(ip), Nquant, quants(ip,:))
+    CALL polygon_properties(dbg, dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),&
+      xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), meanwctrs(ip,:), areas(ip), nvals(ip), xvals(ip),   &
+      mvals(ip), m2vals(ip), stdvals(ip), Nquant, quants(ip,:), nvcoords(ip,:), xvcoords(ip,:),       &
+      meanvnctrs(ip,:), meanvxctrs(ip,:))
   END DO
 
@@ -109 +126 @@
   END SUBROUTINE all_polygons_properties
 
-  SUBROUTINE polygon_properties(dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts, xxtrm,   &
-    yxtrm, meanctr, area, nval, xval, mval, m2val, stdval, Nquant, quant)
+  SUBROUTINE polygon_properties(dbg, dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts,     &
+    xxtrm, yxtrm, meanctr, meanwctr, area, nval, xval, mval, m2val, stdval, Nquant, quant, nvcoord,   &
+    xvcoord, meanvnctr, meanvxctr)
 ! Subroutine to determine the properties of a polygon (as .TRUE. matrix)
 !   Number of grid points
 !   grid-point coordinates of the minimum and maximum of the path along x,y axes
-!   center from the mean of the coordinates of the paths locations
+!   grid coordinates of center from the mean of the coordinates of the poygon locations
+!   lon, lat center from the area weighted mean of the coordinates of the polygon locations
 !   area of the polygon (km2)
 !   minimum and maximum of the values within the polygon
@@ -122 +141 @@
 !   number of quantiles
 !   quantiles of the values within the polygon
+!   grid coordinates of the minimum, maximum value within the polygon
+!   lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
 
   IMPLICIT NONE
 
+  LOGICAL, INTENT(in)                                    :: dbg
   INTEGER, INTENT(in)                                    :: dx, dy, Nquant
   LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: poly
   REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: lon, lat, values
   REAL(r_k), INTENT(in)                                  :: xres, yres
-  CHARACTER(len=50), INTENT(in)                          :: projN
+  CHARACTER(len=1000), INTENT(in)                        :: projN
   INTEGER, INTENT(out)                                   :: Npolypts
   INTEGER, DIMENSION(2), INTENT(out)                     :: xxtrm, yxtrm, meanctr
+  INTEGER, DIMENSION(2), INTENT(out)                     :: nvcoord, xvcoord
+  REAL(r_k), DIMENSION(2), INTENT(out)                   :: meanwctr, meanvnctr, meanvxctr
   REAL(r_k), INTENT(out)                                 :: area
   REAL(r_k), INTENT(out)                                 :: nval, xval, mval, m2val, stdval
@@ -140 +164 @@
   INTEGER                                                :: ierr
   INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: polypts
-  REAL(r_k), DIMENSION(:), ALLOCATABLE                   :: polyvals
+  REAL(r_k), DIMENSION(:), ALLOCATABLE                   :: polyvals, distvn, distvx
+  REAL(r_k)                                              :: lonNADIR, latNADIR
+  REAL(r_k)                                              :: sumRESx, sumRESy
+  REAL(r_k), DIMENSION(dx,dy)                            :: xcorr, ycorr
+  CHARACTER(len=200), DIMENSION(3)                       :: satSvals
+  CHARACTER(len=50)                                      :: projS
+  REAL(r_k)                                              :: sumDISTnlon, sumDISTnlat, sumDISTxlon,   &
+    sumDISTxlat
 
 !!!!!!! Variables
@@ -149 +180 @@
 ! [x/y]res resolution along the x and y axis
 ! projN: name of the projection
+!   'ctsarea': Constant Area
 !   'lon/lat': for regular longitude-latitude
+!   'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR)
 ! Npolypts: number of points of the polygon
 ! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon
 ! meanctr: center from the mean of the coordinates of the polygon
+! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon
 ! area: area of the polygon [km]
 ! [n/x]val: minimum and maximum of the values within the polygon
@@ -159 +193 @@
 ! stdval: standard deviation of the values within the polygon
 ! Nquant: number of quantiles
-! quant: quantiles of the values within the polygon
+! quant: quantiles of the values within the polygon 
+! [n/x]vcoord: grid coordinates of the minimum/maximum of the values within the polygon
+! meanv[n/x]ctr: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon
 
   fname = 'polygon_properties'
+
+  IF (dbg) PRINT *,"  '" // TRIM(fname) // "' ..."
 
   ! Getting grid-point coordinates of the polygon
@@ -176 +214 @@
   CALL ErrMsg(msg, fname, ierr)
 
+  IF (ALLOCATED(distvn)) DEALLOCATE(distvn)
+  ALLOCATE(distvn(Npolypts), STAT=ierr)
+  msg = "Problems allocating 'distvn'"
+  CALL ErrMsg(msg, fname, ierr)
+
+  IF (ALLOCATED(distvx)) DEALLOCATE(distvx)
+  ALLOCATE(distvx(Npolypts), STAT=ierr)
+  msg = "Problems allocating 'distvx'"
+  CALL ErrMsg(msg, fname, ierr)
+
+  IF (projN(1:7) == 'lon/lat') THEN
+    projS = projN
+  ELSE IF (projN(1:7) == 'ctsarea') THEN
+    projS = projN
+  ELSE IF (projN(1:9) == 'nadir-sat') THEN
+    projS = 'nadir-sat'
+    CALL split(projN, ',', 3, satSvals)
+    READ(satSvals(2),'(F200.100)')lonNadir
+    READ(satSvals(3),'(F200.100)')latNadir
+    IF (dbg) PRINT *,"  'nadir-geostationary-satellite' based projection of data with nadir " //      &
+      "location at:", lonNadir, latNadir
+  ELSE
+    msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) // " available ones: " //        &
+       "'ctsarea', 'lon/lat', 'nadir-sat'"
+    CALL ErrMsg(msg,fname,-1)
+  END IF
+
+  area = 0.
+  sumRESx = 0.
+  sumRESy = 0.
+  meanwctr = 0.
+  meanvnctr = 0.
+  meanvxctr = 0.
+
+  nval = fillval64
+  xval = -fillval64
+
   ip = 1
-  area = 0.
   DO i=1,dx
     DO j=1,dy
@@ -183 +257 @@
         polypts(ip,:) = (/ i,j /)
         polyvals(ip) = values(i,j)
-        SELECT CASE (TRIM(projN))
+        SELECT CASE (TRIM(projS))
+          CASE ('ctsarea')
+            ! Constant Area
+            xcorr(i,j) = 1.
+            ycorr(i,j) = 1.
           CASE ('lon/lat')
             ! Area as fixed yres and sinus-lat varying for xres
-            area = area + xres*DABS(DSIN(lon(i,j)*DegRad))*yres
-          CASE DEFAULT
-            msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) //                       &
-              " available ones: 'lon/lat'"
+            xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
+            ycorr(i,j) = 1.
+          CASE ('nadir-sat')
+            ! Area from nadir resolution and degrading as we get far from satellite's nadir
+            ! GOES-E: 0 N, 75 W
+            xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad))
+            ycorr(i,j) = 1.
         END SELECT
+        area = area + xres*xcorr(i,j)*yres*ycorr(i,j)
+        meanwctr(1) = meanwctr(1) + lon(i,j)*xres*xcorr(i,j)
+        meanwctr(2) = meanwctr(2) + lat(i,j)*yres*ycorr(i,j)
+        IF (nval > values(i,j)) THEN
+          nvcoord(1) = i
+          nvcoord(2) = j
+          nval = values(i,j)
+        END IF
+        IF (xval < values(i,j)) THEN
+          xvcoord(1) = i
+          xvcoord(2) = j
+          xval = values(i,j)
+        END IF
         ip = ip + 1
       END IF
     END DO
   END DO
+
+  IF (dbg) THEN
+    PRINT *,'  Values: ', polyvals
+  END IF
+
+  sumRESx = xres*SUM(xcorr)
+  sumRESy = yres*SUM(ycorr)
 
   xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /)
   yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /)
   meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /)
-
+  meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /)
+
+  IF (dbg) THEN
+    PRINT *,'  mean grid center: ', meanctr, ' weighted mean center: ', meanwctr
+  END IF
+  
   ! Statistics of the values within the polygon
   CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
+
+  IF (dbg) THEN
+    PRINT *,'    minimum value: ', nval, ' maximum:', xval, ' mean:', mval
+    PRINT *,'    coor. minimum: ', nvcoord
+    PRINT *,'    coor. maximum: ', xvcoord
+  END IF
+
+  ! Mean center weighted to minimum and maximum values
+  distvn = DABS(polyvals - nval)
+  distvx = DABS(polyvals - xval)
+
+  meanvnctr = 0.
+  meanvxctr = 0.
+  sumDISTnlon = 0.
+  sumDISTnlat = 0.
+  sumDISTxlon = 0.
+  sumDISTxlat = 0.
+
+  ip = 1
+  DO i=1,dx
+    DO j=1,dy
+      IF (poly(i,j)) THEN
+        meanvnctr(1) = meanvnctr(1)+lon(i,j)*distvn(ip)*xres*xcorr(i,j)
+        meanvnctr(2) = meanvnctr(2)+lat(i,j)*distvn(ip)*yres*ycorr(i,j)
+
+        meanvxctr(1) = meanvxctr(1)+lon(i,j)*distvx(ip)*xres*xcorr(i,j)      
+        meanvxctr(2) = meanvxctr(2)+lat(i,j)*distvx(ip)*yres*ycorr(i,j)
+
+        sumDISTnlon = sumDISTnlon + lon(i,j)*distvn(ip)*xres*xcorr(i,j)
+        sumDISTnlat = sumDISTnlat + lat(i,j)*distvn(ip)*yres*ycorr(i,j)
+        sumDISTxlon = sumDISTxlon + lon(i,j)*distvx(ip)*xres*xcorr(i,j)
+        sumDISTxlat = sumDISTxlat + lat(i,j)*distvx(ip)*yres*ycorr(i,j)
+
+        ip = ip + 1
+      END IF
+    END DO
+  END DO
+
+  meanvnctr = (/ meanvnctr(1)/sumDISTnlon, meanvnctr(2)/sumDISTnlat /)
+  meanvxctr = (/ meanvxctr(1)/sumDISTxlon, meanvxctr(2)/sumDISTxlat /)
+
+  IF (dbg) THEN
+    PRINT *,'  mean center to minimum: ', meanvnctr, ' to maximum: ', meanvxctr
+  END IF
 
   ! Quantiles of the values within the polygon