GetR2float Function

private function GetR2float(grid1, grid2, maskReal, maskInteger) result(r2)

compute R2, square of the Pearson correlation coefficient optional argument: mask: compute RMSE only on assigned mask

References:
Karl Pearson (20 June 1895) "Notes on regression and inheritance in the case of two parents," Proceedings of the Royal Society of London, 58 : 240–242

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: grid1
type(grid_real), intent(in) :: grid2
type(grid_real), intent(in), optional :: maskReal
type(grid_integer), intent(in), optional :: maskInteger

Return Value real(kind=float)


Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: den
real(kind=float), public :: den1
real(kind=float), public :: den2
integer(kind=long), public :: i
integer(kind=long), public :: j
real(kind=float), public :: mean1
real(kind=float), public :: mean2
real(kind=float), public :: num

Source Code

FUNCTION GetR2float &
!
(grid1, grid2, maskReal, maskInteger) &
!
RESULT (r2)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid1  
TYPE (grid_real), INTENT(IN) :: grid2 
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
INTEGER (KIND = long) :: i, j
REAL (KIND = float) :: r2, mean1, mean2, num, den1, den2, den
!---------------------------end of declarations--------------------------------

num = 0.
den = 0.
den1 = 0.
den2 = 0.

!check grid1 and grid2 have the same coordinate reference system
 IF ( .NOT. CRSisEqual(grid1,grid2) ) THEN
      CALL Catch ('error', 'GridStatistics',  &
      'calculate R2: ', argument = &
      'coordinate reference system of grid1 differs from grid2' )
END IF

IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate R2: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean of grid1 and grid2
     mean1 = GetMean (grid1, maskReal = maskReal)
     mean2 = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                num = num + ( grid1 % mat (i,j) - mean1 ) * &
                            ( grid2 % mat (i,j) - mean2 )
                den1 = den1 + ( grid1 % mat (i,j) - mean1 ) ** 2.
                den2 = den2 + ( grid2 % mat (i,j) - mean2 ) ** 2.
            END IF
        END DO
    END DO
    
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate R2: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

     !get mean of grid1 and grid2
      mean1 = GetMean (grid1, maskReal = maskReal)
      mean2 = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 num = num + ( grid1 % mat (i,j) - mean1 ) * &
                            ( grid2 % mat (i,j) - mean2 )
                 den1 = den1 + ( grid1 % mat (i,j) - mean1 ) ** 2.
                 den2 = den2 + ( grid2 % mat (i,j) - mean2 ) ** 2.
            END IF
        END DO
    END DO

ELSE !compute over entire grid
   
     !get mean of grid1 and grid2
      mean1 = GetMean (grid1, maskReal = maskReal)
      mean2 = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, grid1 % jdim
        DO i = 1, grid1 % idim
            IF (grid1 % mat(i,j) /= grid1 % nodata) THEN
                 num = num + ( grid1 % mat (i,j) - mean1 ) * &
                            ( grid2 % mat (i,j) - mean2 )
                 den1 = den1 + ( grid1 % mat (i,j) - mean1 ) ** 2.
                 den2 = den2 + ( grid2 % mat (i,j) - mean2 ) ** 2.
            END IF
        END DO
    END DO
END IF

den = (den1 * den2) ** 0.5
IF (den /= 0.) THEN
   r2 = ( num / den ) ** 2.
ELSE
   r2 = -9999.9
END IF


RETURN

END FUNCTION GetR2float