Changeset 1618 in lmdz_wrf


Ignore:
Timestamp:
Sep 5, 2017, 2:52:41 AM (7 years ago)
Author:
lfita
Message:

Adding calculation of properties of the polygons

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/module_scientific.f90

    r1616 r1618  
    33 
    44!!!!!!! Functions & subroutines
     5! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field:
    56! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
    67! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole
    … …  
    1213! paths_border: Subroutine to search the paths of a border field.
    1314! path_properties: Subroutine to determine the properties of a path
     15! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix)
    1416! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
    1517! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
    … …  
    1820! rand_sample: Subroutine to randomly sample a range of indices
    1921! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
     22! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
     23!   series of r_k numbers
    2024! SwapR_K*: Subroutine swaps the values of its two formal arguments.
    2125
    … …  
    2731
    2832  CONTAINS
     33
     34  SUBROUTINE all_polygons_properties(dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN,       &
     35    Npolyptss, xxtrms, yxtrms, meanctrs, areas, nvals, xvals, mvals, m2vals, stdvals, Nquant, quants)
     36! Subroutine to determine the properties of all polygons in a 2D field:
     37!   Number of grid points
     38!   grid-point coordinates of the minimum and maximum of the path along x,y axes
     39!   center from the mean of the coordinates of the paths locations
     40!   area of the polygon (km2)
     41!   minimum and maximum of the values within the polygon
     42!   mean of the values within the polygon
     43!   quadratic mean of the values within the polygon
     44!   standard deviation of the values within the polygon
     45!   number of quantiles
     46!   quantiles of the values within the polygon
     47
     48  IMPLICIT NONE
     49
     50  INTEGER, INTENT(in)                                    :: dx, dy, Npoly, Nquant
     51  INTEGER, DIMENSION(dx,dy), INTENT(in)                  :: polys
     52  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: lon, lat, values
     53  REAL(r_k), INTENT(in)                                  :: xres, yres
     54  CHARACTER(len=50), INTENT(in)                          :: projN
     55  INTEGER, DIMENSION(Npoly), INTENT(out)                 :: Npolyptss
     56  INTEGER, DIMENSION(Npoly, 2), INTENT(out)              :: xxtrms, yxtrms, meanctrs
     57  REAL(r_k), DIMENSION(Npoly), INTENT(out)               :: areas
     58  REAL(r_k), DIMENSION(Npoly), INTENT(out)               :: nvals, xvals, mvals, m2vals, stdvals
     59  REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out)       :: quants
     60
     61! Local
     62  INTEGER                                                :: ip
     63  LOGICAL, DIMENSION(dx,dy)                              :: boolpoly
     64
     65!!!!!!! Variables
     66! dx,dy: size of the space
     67! Npoly: number of polygons
     68! polys: polygon matrix with all polygons (as integer number per polygon)
     69! lon, lat: geographical coordinates of the grid-points of the matrix
     70! values: values of the 2D field to use
     71! [x/y]res resolution along the x and y axis
     72! projN: name of the projection
     73!   'lon/lat': for regular longitude-latitude
     74! Npolyptss: number of points of the polygons
     75! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons
     76! meanctrs: center from the mean of the coordinates of the polygons
     77! areas: area of the polygons [km]
     78! [n/x]vals: minimum and maximum of the values within the polygons
     79! mvals: mean of the values within the polygons
     80! m2vals: quadratic mean of the values within the polygons
     81! stdvals: standard deviation of the values within the polygons
     82! Nquant: number of quantiles
     83! quants: quantiles of the values within the polygons
     84
     85  fname = 'all_polygons_properties'
     86
     87  ! Initializing matrices
     88  Npolyptss = -1
     89  xxtrms = fillval64
     90  yxtrms = fillval64
     91  meanctrs = fillval64
     92  areas = fillval64
     93  nvals = -1
     94  xvals = fillval64
     95  mvals = fillval64
     96  m2vals = fillval64
     97  stdvals = fillval64
     98  quants = fillval64
     99
     100  DO ip=1, Npoly
     101    boolpoly = polys == ip
     102    CALL polygon_properties(dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),     &
     103      xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), areas(ip), nvals(ip), xvals(ip), mvals(ip),         &
     104      m2vals(ip), stdvals(ip), Nquant, quants(ip,:))
     105  END DO
     106
     107  RETURN
     108
     109  END SUBROUTINE all_polygons_properties
     110
     111  SUBROUTINE polygon_properties(dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts, xxtrm,   &
     112    yxtrm, meanctr, area, nval, xval, mval, m2val, stdval, Nquant, quant)
     113! Subroutine to determine the properties of a polygon (as .TRUE. matrix)
     114!   Number of grid points
     115!   grid-point coordinates of the minimum and maximum of the path along x,y axes
     116!   center from the mean of the coordinates of the paths locations
     117!   area of the polygon (km2)
     118!   minimum and maximum of the values within the polygon
     119!   mean of the values within the polygon
     120!   quadratic mean of the values within the polygon
     121!   standard deviation of the values within the polygon
     122!   number of quantiles
     123!   quantiles of the values within the polygon
     124
     125  IMPLICIT NONE
     126
     127  INTEGER, INTENT(in)                                    :: dx, dy, Nquant
     128  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: poly
     129  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: lon, lat, values
     130  REAL(r_k), INTENT(in)                                  :: xres, yres
     131  CHARACTER(len=50), INTENT(in)                          :: projN
     132  INTEGER, INTENT(out)                                   :: Npolypts
     133  INTEGER, DIMENSION(2), INTENT(out)                     :: xxtrm, yxtrm, meanctr
     134  REAL(r_k), INTENT(out)                                 :: area
     135  REAL(r_k), INTENT(out)                                 :: nval, xval, mval, m2val, stdval
     136  REAL(r_k), DIMENSION(Nquant), INTENT(out)              :: quant
     137
     138! Local
     139  INTEGER                                                :: i, j, ip
     140  INTEGER                                                :: ierr
     141  INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: polypts
     142  REAL(r_k), DIMENSION(:), ALLOCATABLE                   :: polyvals
     143
     144!!!!!!! Variables
     145! dx,dy: size of the space
     146! poly: polygon matrix (boolean)
     147! lon, lat: geographical coordinates of the grid-points of the matrix
     148! values: values of the 2D field to use
     149! [x/y]res resolution along the x and y axis
     150! projN: name of the projection
     151!   'lon/lat': for regular longitude-latitude
     152! Npolypts: number of points of the polygon
     153! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon
     154! meanctr: center from the mean of the coordinates of the polygon
     155! area: area of the polygon [km]
     156! [n/x]val: minimum and maximum of the values within the polygon
     157! mval: mean of the values within the polygon
     158! m2val: quadratic mean of the values within the polygon
     159! stdval: standard deviation of the values within the polygon
     160! Nquant: number of quantiles
     161! quant: quantiles of the values within the polygon
     162
     163  fname = 'polygon_properties'
     164
     165  ! Getting grid-point coordinates of the polygon
     166  Npolypts = COUNT(poly)
     167
     168  IF (ALLOCATED(polypts)) DEALLOCATE(polypts)
     169  ALLOCATE(polypts(Npolypts,2), STAT=ierr)
     170  msg = "Problems allocating 'polypts'"
     171  CALL ErrMsg(msg, fname, ierr)
     172
     173  IF (ALLOCATED(polyvals)) DEALLOCATE(polyvals)
     174  ALLOCATE(polyvals(Npolypts), STAT=ierr)
     175  msg = "Problems allocating 'polyvals'"
     176  CALL ErrMsg(msg, fname, ierr)
     177
     178  ip = 1
     179  area = 0.
     180  DO i=1,dx
     181    DO j=1,dy
     182      IF (poly(i,j)) THEN
     183        polypts(ip,:) = (/ i,j /)
     184        polyvals(ip) = values(i,j)
     185        SELECT CASE (TRIM(projN))
     186          CASE ('lon/lat')
     187            ! Area as fixed yres and sinus-lat varying for xres
     188            area = area + xres*DABS(DSIN(lon(i,j)*DegRad))*yres
     189          CASE DEFAULT
     190            msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) //                       &
     191              " available ones: 'lon/lat'"
     192        END SELECT
     193        ip = ip + 1
     194      END IF
     195    END DO
     196  END DO
     197
     198  xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /)
     199  yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /)
     200  meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /)
     201
     202  ! Statistics of the values within the polygon
     203  CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval)
     204
     205  ! Quantiles of the values within the polygon
     206  quant = -9999.d0
     207  IF (Npolypts > Nquant) THEN
     208    CALL quantilesR_K(Npolypts, polyvals, Nquant, quant)
     209  ELSE
     210    CALL SortR_K(polyvals, Npolypts)
     211  END IF
     212
     213  DEALLOCATE (polypts)
     214  DEALLOCATE (polyvals)
     215
     216  RETURN
     217
     218  END SUBROUTINE polygon_properties
    29219
    30220SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
    … …  
    107297  polys = -1
    108298
     299  ! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero)
    109300  Nppt = dx*dy/10
    110301
    … …  
    10911282END SUBROUTINE quantilesR_K
    10921283
     1284
     1285SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev)
     1286! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a
     1287!   series of r_k numbers
     1288
     1289  IMPLICIT NONE
     1290
     1291  INTEGER, INTENT(in)                                    :: Nvals
     1292  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
     1293  REAL(r_k), INTENT(out)                                 :: minv, maxv, mean, mean2, stdev
     1294
     1295!!!!!!! Variables
     1296! Nvals: number of values
     1297! vals: values
     1298! minv: minimum value of values
     1299! maxv: maximum value of values
     1300! mean: mean value of values
     1301! mean2: quadratic mean value of values
     1302! stdev: standard deviation of values
     1303
     1304  minv = MINVAL(vals)
     1305  maxv = MAXVAL(vals)
     1306
     1307  mean=SUM(vals)
     1308  mean2=SUM(vals*vals)
     1309
     1310  mean=mean/Nvals
     1311  mean2=mean2/Nvals
     1312
     1313  stdev=SQRT(mean2 - mean*mean)
     1314
     1315  RETURN
     1316
     1317END SUBROUTINE StatsR_k
     1318
    10931319END MODULE module_scientific
Note: See TracChangeset for help on using the changeset viewer.