Changeset 502 in lmdz_wrf for trunk/tools/validation_sim.py


Ignore:
Timestamp:
Jun 15, 2015, 6:26:59 PM (9 years ago)
Author:
lfita
Message:

Adding Exact time statistics
Adding complementary statistics:
# From:
#Willmott, C. J. 1981. On the validation of models. Physical Geography, 2, 184–194
#Willmott, C. J. (1984). On the evaluation of model performance in physical geography. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds., 443-460
#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R. Legates, J. O'Donnell, and C. M. Rowe (1985), Statistics for the Evaluation and Comparison of Models, J. Geophys. Res., 90(C5), 8995-9005
#Legates, D. R., and G. J. McCabe? Jr. (1999), Evaluating the Use of "Goodness-of-Fit" Measures in Hydrologic and Hydroclimatic Model Validation, Water Resour. Res., 35(1), 233–241
#
# Deviation of residuals (SDR)
# Index of Efficiency (IOE)
# Index of Agreement (IOA)
# Fractional Mean Bias (FMB)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/validation_sim.py

    r501 r502  
    12681268(opts, args) = parser.parse_args()
    12691269
     1270####### ###### ##### #### ### ## #
     1271# Number of different statistics according to the temporal coincidence
     1272#  0: Exact time
     1273#  1: Simulation values between consecutive observed times
     1274Nstsim = 2
     1275
     1276stdescsim = ['E', 'B']
     1277Lstdescsim = ['exact time', 'between observational time-steps']
     1278
    12701279#######    #######
    12711280## MAIN
    … …  
    15941603newdim = onewnc.createDimension('stats',5)
    15951604newdim = onewnc.createDimension('tstats',6)
     1605newdim = onewnc.createDimension('Nstsim, 2)
    15961606
    15971607# Variable dimensions
    … …  
    16181628basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
    16191629writing_str_nc(newvar, ostatsn, StringLength)
     1630
     1631newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
     1632basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
     1633writing_str_nc(newvar, Lstdescsim, StringLength)
     1634
    16201635
    16211636if obskind == 'trajectory':
    … …  
    16711686        oFillValue = None
    16721687
    1673 # Observed values temporally around coincident times
    1674 #    simobsvalues, simobsSvalues, simobsTvalues, simobsTtvalues, trjsim =             \
    1675 #        getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, ovobs,      \
    1676 #        exacttvalues, oFillValue, Ngrid)
     1688# Observed values temporally exact times
     1689    Esimobsvalues, EsimobsSvalues, EsimobsTvalues, EsimobsTtvalues, trjsim =         \
     1690        getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, ovobs,      \
     1691        exacttvalues, oFillValue, Ngrid)
    16771692
    16781693# Observed values temporally around coincident times
    … …  
    16821697
    16831698# Re-arranging values...
     1699    Earrayvals = np.array(Esimobsvalues)
    16841700    arrayvals = np.array(simobsvalues)
    16851701    if len(valvars[ivar]) > 2:
    16861702        const=np.float(valvars[ivar][3])
    16871703        if valvars[ivar][2] == 'sumc':
     1704            EsimobsSvalues = EsimobsSvalues + const
     1705            Earrayvals[:,0] = Earrayvals[:,0] + const
    16881706            simobsSvalues = simobsSvalues + const
    16891707            arrayvals[:,0] = arrayvals[:,0] + const
    16901708        elif valvars[ivar][2] == 'subc':
     1709            EsimobsSvalues = EsimobsSvalues - const
     1710            Earrayvals[:,0] = Earrayvals[:,0] - const
    16911711            simobsSvalues = simobsSvalues - const
    16921712            arrayvals[:,0] = arrayvals[:,0] - const
    16931713        elif valvars[ivar][2] == 'mulc':
     1714            EsimobsSvalues = EsimobsSvalues * const
     1715            Earrayvals[:,0] = Earrayvals[:,0] * const
    16941716            simobsSvalues = simobsSvalues * const
    16951717            arrayvals[:,0] = arrayvals[:,0] * const
    16961718        elif valvars[ivar][2] == 'divc':
     1719            EsimobsSvalues = EsimobsSvalues / const
     1720            Earrayvals[:,0] = Earrayvals[:,0] / const
    16971721            simobsSvalues = simobsSvalues / const
    16981722            arrayvals[:,0] = arrayvals[:,0] / const
    … …  
    17031727
    17041728# statisics sim
    1705     simstats = np.zeros((5), dtype=np.float)
    1706     simstats[0] = np.min(arrayvals[:,0])
    1707     simstats[1] = np.max(arrayvals[:,0])
    1708     simstats[2] = np.mean(arrayvals[:,0])
    1709     simstats[3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
    1710     simstats[4] = np.sqrt(simstats[3] - simstats[2]*simstats[2])
     1729    simstats = np.zeros((2,5), dtype=np.float)
     1730    simstats[0,0] = np.min(Earrayvals[:,0])
     1731    simstats[0,1] = np.max(Earrayvals[:,0])
     1732    simstats[0,2] = np.mean(Earrayvals[:,0])
     1733    simstats[0,3] = np.mean(Earrayvals[:,0]*Earrayvals[:,0])
     1734    simstats[0,4] = np.sqrt(simstats[0,3] - simstats[0,2]*simstats[0,2])
     1735
     1736    simstats[1,0] = np.min(arrayvals[:,0])
     1737    simstats[1,1] = np.max(arrayvals[:,0])
     1738    simstats[1,2] = np.mean(arrayvals[:,0])
     1739    simstats[1,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
     1740    simstats[1,4] = np.sqrt(simstats[1,3] - simstats[1,2]*simstats[1,2])
    17111741
    17121742# statisics obs
    … …  
    17241754 
    17251755# Statistics sim-obs
    1726     simobsstats = np.zeros((9), dtype=np.float)
     1756    simobsstats = np.zeros((2,13), dtype=np.float)
     1757    Ediffvals = np.zeros((Nexactt), dtype=np.float)
    17271758    diffvals = np.zeros((Ncoindt), dtype=np.float)
    17281759
     1760    Ediffvals = Earrayvals[:,0] - obsmask
    17291761    diffvals = arrayvals[:,0] - obsmask
    17301762
    1731     simobsstats[0] = simstats[0] - obsstats[0]
    1732     simobsstats[1] = np.mean(arrayvals[:,0]*obsmask)
    1733     simobsstats[2] = np.min(diffvals)
    1734     simobsstats[3] = np.max(diffvals)
    1735     simobsstats[4] = np.mean(diffvals)
    1736     simobsstats[5] = np.mean(np.abs(diffvals))
    1737     simobsstats[6] = np.sqrt(np.mean(diffvals*diffvals))
    1738     simobsstats[7], simobsstats[8] = sts.mstats.pearsonr(arrayvals[:,0], arrayvals[:,1])
     1763    Ediff2vals = Ediffvals*Ediffvals
     1764    Esumdiff = Ediffvals.sum()
     1765    Esumdiff2 = Ediff2vals.sum()
     1766
     1767    diff2vals = diffvals*diffvals
     1768    sumdiff = diffvals.sum()
     1769    sumdiff2 = diff2vals.sum()
     1770
     1771    simobsstats[0,0] = simstats[0,0] - obsstats[0]
     1772    simobsstats[0,1] = np.mean(Earrayvals[:,0]*obsmask)
     1773    simobsstats[0,2] = Ediffvals.min()
     1774    simobsstats[0,3] = Ediffvals.max()
     1775    simobsstats[0,4] = Ediffvals.mean()
     1776    simobsstats[0,5] = np.abs(Ediffvals).mean()
     1777    simobsstats[0,6] = np.sqrt(Ediff2vals.mean())
     1778    simobsstats[0,7], simobsstats[0,8] = sts.mstats.pearsonr(Earrayvals[:,0],        \
     1779      obsmask)
     1780
     1781    simobsstats[1,0] = simstats[1,0] - obsstats[0]
     1782    simobsstats[1,1] = np.mean(arrayvals[:,0]*obsmask)
     1783    simobsstats[1,2] = diffvals.min()
     1784    simobsstats[1,3] = diffvals.max()
     1785    simobsstats[1,4] = diffvals.mean()
     1786    simobsstats[1,5] = np.abs(diffvals).mean()
     1787    simobsstats[1,6] = np.sqrt(diff2vals.mean())
     1788    simobsstats[1,7], simobsstats[1,8] = sts.mstats.pearsonr(arrayvals[:,0],         \
     1789      obsmask)
     1790# From:
     1791#Willmott, C. J. 1981. On the validation of models. Physical Geography, 2, 184–194
     1792#Willmott, C. J. (1984). On the evaluation of model performance in physical geography. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds., 443-460
     1793#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R. Legates, J. O'Donnell, and C. M. Rowe (1985), Statistics for the Evaluation and Comparison of Models, J. Geophys. Res., 90(C5), 8995-9005
     1794#Legates, D. R., and G. J. McCabe Jr. (1999), Evaluating the Use of "Goodness-of-Fit" Measures in Hydrologic and Hydroclimatic Model Validation, Water Resour. Res., 35(1), 233–241
     1795#
     1796# Deviation of residuals (SDR)
     1797    simobsstats[0,9] = np.mean(np.sqrt(np.abs((Ediffvals-simobsstats[0,0])*(Ediffvals-\
     1798      obsstats[0]))))
     1799    simobsstats[1,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[1,0])*(diffvals- \
     1800      obsstats[0]))))
     1801# Index of Efficiency (IOE)
     1802    obsbias2series = (obsmask - obsstats[0])*(obsmask - obsstats[0])
     1803    sumobsbias2series = obsbias2series.sum()
     1804
     1805    simobsstats[1,10] = 1. - Esumdiff2/(sumobsbias2series)
     1806    simobsstats[1,10] = 1. - sumdiff2/(sumobsbias2series)
     1807# Index of Agreement (IOA)
     1808    Esimbias2series = Earrayvals[:,0] - obsstats[0]
     1809    simbias2series = arrayvals[:,0] - obsstats[0]
     1810
     1811    obsbias2series = obsmask - obsstats[0]
     1812
     1813    Eabssimbias2series = np.abs(Esimbias2series)
     1814    abssimbias2series = np.abs(simbias2series)
     1815    absobsbias2series = obbias2series.abs()
     1816
     1817    Eabssimobsbias2series = (Eabssimbias2series+absobsbias2series)*(Eabssimbias2series +\
     1818      absobsbias2series)
     1819    abssimobsbias2series = (abssimbias2series+absobsbias2series)*(abssimbias2series +\
     1820      absobsbias2series)
     1821
     1822    simobsstats[0,11] = 1. - Esumdiff2/(Eabssimobsbias2series.sum())
     1823    simobsstats[1,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
     1824# Fractional Mean Bias (FMB)
     1825    simobsstats[0,12] = (simstats[0,0]-obsstats[0])/(0.5*(simstats[0,0]+obsstats[0]))
     1826    simobsstats[1,12] = (simstats[1,0]-obsstats[0])/(0.5*(simstats[1,0]+obsstats[0]))
    17391827
    17401828# Statistics around sim values
    … …  
    17931881
    17941882# sim Statistics
     1883    for ist in range(Nstsim):
    17951884    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
    1796         newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('stats'),   \
    1797           fill_value=fillValueF)
     1885        newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('Nstsim',   \
     1886          'stats'), fill_value=fillValueF)
    17981887        descvar = 'simulated statisitcs: ' + valvars[ivar][0]
    17991888        basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, ovobs.getncattr('units'))
    … …  
    18111900# sim-obs Statistics
    18121901    if not searchInlist(onewnc.variables,varsimobs + 'st'):
    1813         newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('gstats'),            \
     1902        newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('Nstsim', 'gstats'),  \
    18141903          fill_value=fillValueF)
    18151904        descvar = 'simulated-observed statisitcs: ' + varsimobs
Note: See TracChangeset for help on using the changeset viewer.