PSM for tree TRW (VS-Lite)

In this tutorial, we introduce the PSM for tree TRW (VS-Lite) in cfr.

%load_ext autoreload
%autoreload 2

import cfr
print(cfr.__version__)

Data preparation

Proxy

pdb = cfr.ProxyDatabase().fetch('PAGES2kv2')

pobj = pdb.records['NAm_018']
fig, ax = pobj.plot()

Model

model_tas = cfr.ClimateField().fetch('iCESM_past1000historical/tas')
model_pr = cfr.ClimateField().fetch('iCESM_past1000historical/pr')

>>> The target file seems existed at: ./data/tas_sfc_Amon_iCESM_past1000historical_085001-200512.nc . Loading from it instead of downloading ...
>>> The target file seems existed at: ./data/pr_sfc_Amon_iCESM_past1000historical_085001-200512.nc . Loading from it instead of downloading ...

model_tas.da

<xarray.DataArray 'tas' (time: 13872, lat: 96, lon: 144)>
[191766528 values with dtype=float32]
Coordinates:
  * time     (time) object 0850-01-17 00:00:00 ... 2005-12-17 00:00:00
  * lat      (lat) float32 -90.0 -88.11 -86.21 -84.32 ... 84.32 86.21 88.11 90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    long_name:  Reference height temperature
    units:      K

xarray.DataArray

'tas'

• time: 13872
• lat: 96
• lon: 144

• ...

  [191766528 values with dtype=float32]

• Coordinates: (3)
  • time
    (time)
    object
    0850-01-17 00:00:00 ... 2005-12-...

    description :

    time

    long_name :

    time

    standard_name :

    time

    array([cftime.DatetimeNoLeap(850, 1, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(850, 2, 14, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(850, 3, 17, 0, 0, 0, 0, has_year_zero=True), ...,
           cftime.DatetimeNoLeap(2005, 10, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2005, 11, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2005, 12, 17, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • lat
    (lat)
    float32
    -90.0 -88.11 -86.21 ... 88.11 90.0

    description :

    latitude

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    array([-90.      , -88.10526 , -86.210526, -84.31579 , -82.42105 , -80.52631 ,
           -78.63158 , -76.73684 , -74.8421  , -72.947365, -71.052635, -69.1579  ,
           -67.26316 , -65.36842 , -63.473682, -61.57895 , -59.68421 , -57.789474,
           -55.894737, -54.      , -52.105263, -50.210526, -48.31579 , -46.42105 ,
           -44.526318, -42.63158 , -40.736843, -38.842106, -36.94737 , -35.05263 ,
           -33.157894, -31.263159, -29.368422, -27.473684, -25.578947, -23.68421 ,
           -21.789474, -19.894737, -18.      , -16.105263, -14.210526, -12.315789,
           -10.421053,  -8.526316,  -6.631579,  -4.736842,  -2.842105,  -0.947368,
             0.947368,   2.842105,   4.736842,   6.631579,   8.526316,  10.421053,
            12.315789,  14.210526,  16.105263,  18.      ,  19.894737,  21.789474,
            23.68421 ,  25.578947,  27.473684,  29.368422,  31.263159,  33.157894,
            35.05263 ,  36.94737 ,  38.842106,  40.736843,  42.63158 ,  44.526318,
            46.42105 ,  48.31579 ,  50.210526,  52.105263,  54.      ,  55.894737,
            57.789474,  59.68421 ,  61.57895 ,  63.473682,  65.36842 ,  67.26316 ,
            69.1579  ,  71.052635,  72.947365,  74.8421  ,  76.73684 ,  78.63158 ,
            80.52631 ,  82.42105 ,  84.31579 ,  86.210526,  88.10526 ,  90.      ],
          dtype=float32)

  • lon
    (lon)
    float32
    0.0 2.5 5.0 ... 352.5 355.0 357.5

    description :

    longitude

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
            25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
            50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
            75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
           100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
           125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
           150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
           175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
           200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
           275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
           300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
           325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
           350. , 352.5, 355. , 357.5], dtype=float32)

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(CFTimeIndex([0850-01-17 00:00:00, 0850-02-14 00:00:00, 0850-03-17 00:00:00,
                 0850-04-16 00:00:00, 0850-05-17 00:00:00, 0850-06-16 00:00:00,
                 0850-07-17 00:00:00, 0850-08-17 00:00:00, 0850-09-16 00:00:00,
                 0850-10-17 00:00:00,
                 ...
                 2005-03-17 00:00:00, 2005-04-16 00:00:00, 2005-05-17 00:00:00,
                 2005-06-16 00:00:00, 2005-07-17 00:00:00, 2005-08-17 00:00:00,
                 2005-09-16 00:00:00, 2005-10-17 00:00:00, 2005-11-16 00:00:00,
                 2005-12-17 00:00:00],
                dtype='object', length=13872, calendar='noleap', freq='None'))

  • lat
    PandasIndex

    PandasIndex(Float64Index([              -90.0,  -88.10526275634766,  -86.21052551269531,
                   -84.31578826904297,  -82.42105102539062,  -80.52631378173828,
                   -78.63157653808594,   -76.7368392944336,  -74.84210205078125,
                    -72.9473648071289,   -71.0526351928711,  -69.15789794921875,
                    -67.2631607055664,  -65.36842346191406,  -63.47368240356445,
                  -61.578948974609375,  -59.68421173095703,  -57.78947448730469,
                  -55.894737243652344,               -54.0, -52.105262756347656,
                   -50.21052551269531,  -48.31578826904297, -46.421051025390625,
                   -44.52631759643555,   -42.6315803527832,  -40.73684310913086,
                  -38.842105865478516,  -36.94736862182617,  -35.05263137817383,
                  -33.157894134521484, -31.263158798217773,  -29.36842155456543,
                  -27.473684310913086, -25.578947067260742,   -23.6842098236084,
                  -21.789474487304688, -19.894737243652344,               -18.0,
                  -16.105262756347656, -14.210526466369629, -12.315789222717285,
                  -10.421052932739258,  -8.526315689086914, -6.6315789222717285,
                   -4.736842155456543, -2.8421051502227783, -0.9473684430122375,
                   0.9473684430122375,  2.8421051502227783,   4.736842155456543,
                   6.6315789222717285,   8.526315689086914,  10.421052932739258,
                   12.315789222717285,  14.210526466369629,  16.105262756347656,
                                 18.0,  19.894737243652344,  21.789474487304688,
                     23.6842098236084,  25.578947067260742,  27.473684310913086,
                    29.36842155456543,  31.263158798217773,  33.157894134521484,
                    35.05263137817383,   36.94736862182617,  38.842105865478516,
                    40.73684310913086,    42.6315803527832,   44.52631759643555,
                   46.421051025390625,   48.31578826904297,   50.21052551269531,
                   52.105262756347656,                54.0,  55.894737243652344,
                    57.78947448730469,   59.68421173095703,  61.578948974609375,
                    63.47368240356445,   65.36842346191406,    67.2631607055664,
                    69.15789794921875,    71.0526351928711,    72.9473648071289,
                    74.84210205078125,    76.7368392944336,   78.63157653808594,
                    80.52631378173828,   82.42105102539062,   84.31578826904297,
                    86.21052551269531,   88.10526275634766,                90.0],
                 dtype='float64', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Float64Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,
                   22.5,
                  ...
                  335.0, 337.5, 340.0, 342.5, 345.0, 347.5, 350.0, 352.5, 355.0,
                  357.5],
                 dtype='float64', name='lon', length=144))

• Attributes: (2)

  long_name :

  Reference height temperature

  units :

  K

Instrumental observation

obs_tas = cfr.ClimateField().fetch('CRUTSv4.07/tas', vn='tmp')
obs_pr = cfr.ClimateField().fetch('CRUTSv4.07/pr', vn='pre')

>>> The target file seems existed at: ./data/cru_ts4.07.1901.2022.tmp.dat.nc.gz . Loading from it instead of downloading ...
>>> The target file seems existed at: ./data/cru_ts4.07.1901.2022.pre.dat.nc.gz . Loading from it instead of downloading ...

obs_pr = obs_pr.rename('pr')
obs_tas = obs_tas.rename('tas')

Get climate data for a specific ProxyRecord

%%time

pobj.del_clim()
pobj.get_clim(model_tas, tag='model')
pobj.get_clim(model_pr, tag='model')
pobj.get_clim(obs_tas, tag='obs')
pobj.get_clim(obs_pr, tag='obs')

CPU times: user 8.01 ms, sys: 152 ms, total: 160 ms
Wall time: 1.23 s

pobj.clim['obs.tas'].da

<xarray.DataArray 'tas' (time: 1464)>
array([ 0.8       ,  0.90000004,  2.2       , ..., 11.1       ,
        1.7       ,  0.3       ], dtype=float32)
Coordinates:
    lon      float32 241.8
    lat      float32 36.25
  * time     (time) object 1901-01-16 00:00:00 ... 2022-12-16 00:00:00
Attributes:
    long_name:                   near-surface temperature
    units:                       degrees Celsius
    correlation_decay_distance:  1200.0

xarray.DataArray

'tas'

• time: 1464

• 0.8 0.9 2.2 3.4 8.4 14.0 19.3 ... 16.8 20.4 20.4 17.8 11.1 1.7 0.3

  array([ 0.8       ,  0.90000004,  2.2       , ..., 11.1       ,
          1.7       ,  0.3       ], dtype=float32)

• Coordinates: (3)
  • lon
    ()
    float32
    241.8

    long_name :

    longitude

    units :

    degrees_east

    array(241.75, dtype=float32)

  • lat
    ()
    float32
    36.25

    long_name :

    latitude

    units :

    degrees_north

    array(36.25, dtype=float32)

  • time
    (time)
    object
    1901-01-16 00:00:00 ... 2022-12-...

    long_name :

    time

    array([cftime.DatetimeGregorian(1901, 1, 16, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1901, 2, 15, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1901, 3, 16, 0, 0, 0, 0, has_year_zero=False),
           ...,
           cftime.DatetimeGregorian(2022, 10, 16, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2022, 11, 16, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2022, 12, 16, 0, 0, 0, 0, has_year_zero=False)],
          dtype=object)

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1901-01-16 00:00:00, 1901-02-15 00:00:00, 1901-03-16 00:00:00,
                 1901-04-16 00:00:00, 1901-05-16 00:00:00, 1901-06-16 00:00:00,
                 1901-07-16 00:00:00, 1901-08-16 00:00:00, 1901-09-16 00:00:00,
                 1901-10-16 00:00:00,
                 ...
                 2022-03-16 00:00:00, 2022-04-16 00:00:00, 2022-05-16 00:00:00,
                 2022-06-16 00:00:00, 2022-07-16 00:00:00, 2022-08-16 00:00:00,
                 2022-09-16 00:00:00, 2022-10-16 00:00:00, 2022-11-16 00:00:00,
                 2022-12-16 00:00:00],
                dtype='object', length=1464, calendar='standard', freq='None'))

• Attributes: (3)

  long_name :

  near-surface temperature

  units :

  degrees Celsius

  correlation_decay_distance :

  1200.0

Create a PSM object

mdl = cfr.psm.VSLite(pobj)

%%time

mdl.calibrate()  # octave is required. For example, to install on macOS: `brew install octave`
print(mdl.calib_details)

warning: function /Users/fengzhu/Apps/miniconda3/envs/cfr-env/lib/python3.9/site-packages/PyVSL/rng.m shadows a core library function
warning: called from
    _pyeval at line 57 column 30




Working on chain 1 out of 3...
Working on chain 2 out of 3...
Working on chain 3 out of 3...
warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 456 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 457 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 458 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 459 column 9
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 461 column 13
    _pyeval at line 57 column 30

warning: implicit conversion from numeric to char
warning: called from
    estimate_vslite_params_v2_3 at line 466 column 9
    _pyeval at line 57 column 30

    Rhat for T1, T2, M1, M2, sigma2rw:
   1.0000   1.0022   1.0016   1.0004   1.0119
{'T1': 5.661133563105939, 'T2': 18.10182563306548, 'M1': 0.0352692544593921, 'M2': 0.27955590028116684, 'Rhats': array([[1.00001939, 1.00215255, 1.00161344, 1.00035685, 1.01190535]]), 'convwarning': 0.0}
CPU times: user 36 ms, sys: 37.9 ms, total: 73.9 ms
Wall time: 29.5 s

%%time
pp = mdl.forward()

CPU times: user 181 ms, sys: 4.46 ms, total: 185 ms
Wall time: 184 ms

fig, ax = pp.plot()