import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import esmlab
import warnings
# supress Runtime warnings
warnings.simplefilter("ignore", category=RuntimeWarning)
# Turn warmings back on
#warnings.resetwarnings()
# This is a 4x5 I2000 control simulation with GSWP forcing (cycled 2000-2004) for 30 years
scr_dir = '/glade/scratch/kdagon/archive/'
case_name = 'hydro_ensemble_default_params'
file_path = scr_dir+case_name+'/lnd/hist/'
print(file_path)
/glade/scratch/kdagon/archive/hydro_ensemble_default_params/lnd/hist/
model_yr = 30
fils_in = xr.open_mfdataset(file_path+'*.clm2.h0.00'+str(model_yr)+'*.nc', combine='by_coords')
fils_in
<xarray.Dataset> Dimensions: (cft: 2, glc_nec: 10, hist_interval: 2, lat: 46, levdcmp: 1, levgrnd: 25, levlak: 10, levsoi: 20, lon: 72, ltype: 9, natpft: 15, nvegwcs: 4, time: 12) Coordinates: * lon (lon) float32 0.0 5.0 10.0 ... 345.0 350.0 355.0 * lat (lat) float32 -90.0 -86.0 -82.0 ... 82.0 86.0 90.0 * levdcmp (levdcmp) float32 1.0 * levgrnd (levgrnd) float32 0.01 0.04 ... 28.870724 41.998436 * levlak (levlak) float32 0.05 0.6 2.1 ... 25.6 34.325 44.775 * time (time) object 0030-02-01 00:00:00 ... 0031-01-01 00:00:00 Dimensions without coordinates: cft, glc_nec, hist_interval, levsoi, ltype, natpft, nvegwcs Data variables: mcdate (time) int32 dask.array<chunksize=(1,), meta=np.ndarray> mcsec (time) int32 dask.array<chunksize=(1,), meta=np.ndarray> mdcur (time) int32 dask.array<chunksize=(1,), meta=np.ndarray> mscur (time) int32 dask.array<chunksize=(1,), meta=np.ndarray> nstep (time) int32 dask.array<chunksize=(1,), meta=np.ndarray> time_bounds (time, hist_interval) object dask.array<chunksize=(1, 2), meta=np.ndarray> date_written (time) |S16 dask.array<chunksize=(1,), meta=np.ndarray> time_written (time) |S16 dask.array<chunksize=(1,), meta=np.ndarray> area (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> landfrac (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> landmask (time, lat, lon) float64 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> pftmask (time, lat, lon) float64 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> nbedrock (time, lat, lon) float64 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ATM_TOPO (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> BCDEP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> BTRAN2 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> BTRANMN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> DSL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> DSTDEP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> DSTFLXT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> EFLXBUILD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> EFLX_DYNBAL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> EFLX_GRND_LAKE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> EFLX_LH_TOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> EFLX_LH_TOT_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ELAI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ERRH2O (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ERRH2OSNO (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ERRSEB (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ERRSOI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ERRSOL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ESAI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FCEV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FCOV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FCTR (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FGEV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FGR (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FGR12 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FH2OSFC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FIRA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FIRA_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FIRE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FIRE_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FLDS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FPSN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSAT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSND (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSNDLN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSNI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSVD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSVDLN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSVI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSDSVILN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH_G (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH_PRECIP_CONVERSION (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH_RUNOFF_ICE_TO_LIQ (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH_TO_COUPLER (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSH_V (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSM (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSNO (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSNO_EFF (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSR (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSRND (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSRNDLN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSRNI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSRVD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSRVDLN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> FSRVI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> GSSHA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> GSSUN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> H2OCAN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> H2OSFC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> H2OSNO (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> H2OSNO_TOP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> H2OSOI (time, levsoi, lat, lon) float32 dask.array<chunksize=(1, 20, 46, 72), meta=np.ndarray> HEAT_CONTENT1 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HEAT_FROM_AC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HIA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HIA_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HIA_U (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HUMIDEX (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HUMIDEX_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> HUMIDEX_U (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ICE_CONTENT1 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> JMX25T (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> Jmx25Z (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LAISHA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LAISUN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LAKEICEFRAC_SURF (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LAKEICETHICK (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LIQCAN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LIQUID_CONTENT1 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> LNC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_acetaldehyde (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_acetic_acid (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_acetone (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_carene_3 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_ethanol (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_formaldehyde (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_isoprene (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_methanol (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_pinene_a (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> MEG_thujene_a (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> OCDEP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> PARVEGLN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> PBOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> PCO2 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> PCT_CFT (time, cft, lat, lon) float32 dask.array<chunksize=(1, 2, 46, 72), meta=np.ndarray> PCT_GLC_MEC (time, glc_nec, lat, lon) float32 dask.array<chunksize=(1, 10, 46, 72), meta=np.ndarray> PCT_LANDUNIT (time, ltype, lat, lon) float32 dask.array<chunksize=(1, 9, 46, 72), meta=np.ndarray> PCT_NAT_PFT (time, natpft, lat, lon) float32 dask.array<chunksize=(1, 15, 46, 72), meta=np.ndarray> Q2M (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QBOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QCHARGE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QDRAI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QDRAI_PERCH (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QDRAI_XS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QDRIP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLOOD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_DEW_GRND (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_DEW_SNOW (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_EVAP_TOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_ICE_DYNBAL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_LIQ_DYNBAL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_SNOW_DRAIN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_SNOW_DRAIN_ICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QFLX_SUB_SNOW (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QH2OSFC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QICE_FRZ (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QICE_MELT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QINFL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QINTR (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QIRRIG (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QOVER (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QRGWL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QRUNOFF (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QRUNOFF_ICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QRUNOFF_ICE_TO_COUPLER (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QRUNOFF_TO_COUPLER (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNOCPLIQ (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNOEVAP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNOFRZ (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNOFRZ_ICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNOMELT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNOMELT_ICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNO_TEMPUNLOAD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNO_WINDUNLOAD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSNWCPICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSOIL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QSOIL_ICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QVEGE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> QVEGT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> RAIN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> RAIN_FROM_ATM (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> RH2M (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> RSSHA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> RSSUN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SABG (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SABG_PEN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SABV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SMP (time, levgrnd, lat, lon) float32 dask.array<chunksize=(1, 25, 46, 72), meta=np.ndarray> SNOBCMCL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOBCMSL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOCAN (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNODSTMCL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNODSTMSL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOFSRND (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOFSRNI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOFSRVD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOFSRVI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOINTABS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOOCMCL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOOCMSL (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOTXMASS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOUNLOAD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOW (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOWDP (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOWICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOWLIQ (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOW_DEPTH (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOW_FROM_ATM (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOW_PERSISTENCE (time, lat, lon) timedelta64[ns] dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOW_SINKS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SNOW_SOURCES (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SOILICE (time, levsoi, lat, lon) float32 dask.array<chunksize=(1, 20, 46, 72), meta=np.ndarray> SOILLIQ (time, levsoi, lat, lon) float32 dask.array<chunksize=(1, 20, 46, 72), meta=np.ndarray> SOILRESIS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SOILWATER_10CM (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SWBGT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SWBGT_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> SWBGT_U (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TAUX (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TAUY (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TBOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TBUILD (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TG (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TH2OSFC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> THBOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TKE1 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TLAI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TLAKE (time, levlak, lat, lon) float32 dask.array<chunksize=(1, 10, 46, 72), meta=np.ndarray> TOTSOILICE (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TOTSOILLIQ (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TPU25T (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TREFMNAV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TREFMXAV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TSA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TSAI (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TSOI (time, levgrnd, lat, lon) float32 dask.array<chunksize=(1, 25, 46, 72), meta=np.ndarray> TSOI_10CM (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TSOI_ICE (time, levgrnd, lat, lon) float32 dask.array<chunksize=(1, 25, 46, 72), meta=np.ndarray> TV (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> TWS (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> U10 (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> U10_DUST (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> URBAN_AC (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> URBAN_HEAT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> VCMX25T (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> VEGWP (time, nvegwcs, lat, lon) float32 dask.array<chunksize=(1, 4, 46, 72), meta=np.ndarray> VOLR (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> VOLRMCH (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> Vcmx25Z (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> WA (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> WASTEHEAT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> WBT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> WBT_R (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> WBT_U (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> WIND (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ZBOT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ZWT (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> ZWT_PERCH (time, lat, lon) float32 dask.array<chunksize=(1, 46, 72), meta=np.ndarray> Attributes: title: CLM History file information comment: NOTE: None of the variables ar... Conventions: CF-1.0 history: created on 01/24/20 17:02:25 source: Community Land Model CLM4.0 hostname: cheyenne username: kdagon version: unknown revision_id: $Id: histFileMod.F90 42903 201... case_title: UNSET case_id: hydro_ensemble_default_params Surface_dataset: surfdata_4x5_16pfts_Irrig_CMIP... Initial_conditions_dataset: finidat_interp_dest.nc PFT_physiological_constants_dataset: clm5_params.c171117.nc ltype_vegetated_or_bare_soil: 1 ltype_crop: 2 ltype_UNUSED: 3 ltype_landice_multiple_elevation_classes: 4 ltype_deep_lake: 5 ltype_wetland: 6 ltype_urban_tbd: 7 ltype_urban_hd: 8 ltype_urban_md: 9 ctype_vegetated_or_bare_soil: 1 ctype_crop: 2 ctype_crop_noncompete: 2*100+m, m=cft_lb,cft_ub ctype_landice: 3 ctype_landice_multiple_elevation_classes: 4*100+m, m=1,glcnec ctype_deep_lake: 5 ctype_wetland: 6 ctype_urban_roof: 71 ctype_urban_sunwall: 72 ctype_urban_shadewall: 73 ctype_urban_impervious_road: 74 ctype_urban_pervious_road: 75 cft_c3_crop: 1 cft_c3_irrigated: 2 time_period_freq: month_1 Time_constant_3Dvars_filename: ./hydro_ensemble_default_param... Time_constant_3Dvars: ZSOI:DZSOI:WATSAT:SUCSAT:BSW:H...
TSA = fils_in.TSA # this is now a DataArray of 2m air temperature
TSA
<xarray.DataArray 'TSA' (time: 12, lat: 46, lon: 72)> dask.array<concatenate, shape=(12, 46, 72), dtype=float32, chunksize=(1, 46, 72), chunktype=numpy.ndarray> Coordinates: * lon (lon) float32 0.0 5.0 10.0 15.0 20.0 ... 340.0 345.0 350.0 355.0 * lat (lat) float32 -90.0 -86.0 -82.0 -78.0 -74.0 ... 78.0 82.0 86.0 90.0 * time (time) object 0030-02-01 00:00:00 ... 0031-01-01 00:00:00 Attributes: long_name: 2m air temperature units: K cell_methods: time: mean
TSA_amean = TSA.mean(dim='time')
TSA_amean.plot(robust=True);
# Just pull the first time index since these are time-invariant fields
area = fils_in.area[0]
landfrac = fils_in.landfrac[0]
arwt = area * landfrac
# re-assign long_name attribute
arwt.attrs['long_name'] = 'grid cell areas multiplied by landfrac'
# quick plot of area weights
arwt.plot(robust=True);
# normalize by sum of area weights
arwt_norm_sum = arwt/arwt.sum()
arwt_norm_sum.attrs['long_name'] = 'area weights normalized by sum'
arwt_norm_sum.attrs['units'] = ''
arwt_norm_sum
<xarray.DataArray (lat: 46, lon: 72)> dask.array<truediv, shape=(46, 72), dtype=float32, chunksize=(46, 72), chunktype=numpy.ndarray> Coordinates: * lon (lon) float32 0.0 5.0 10.0 15.0 20.0 ... 340.0 345.0 350.0 355.0 * lat (lat) float32 -90.0 -86.0 -82.0 -78.0 -74.0 ... 78.0 82.0 86.0 90.0 time object 0030-02-01 00:00:00 Attributes: long_name: area weights normalized by sum units:
arwt_norm_sum.plot(robust=True);
# normalize by maximum of area weights
arwt_norm_max = arwt/arwt.max()
arwt_norm_max.attrs['long_name'] = 'area weights normalized by max'
arwt_norm_max.attrs['units'] = ''
arwt_norm_max
<xarray.DataArray (lat: 46, lon: 72)> dask.array<truediv, shape=(46, 72), dtype=float32, chunksize=(46, 72), chunktype=numpy.ndarray> Coordinates: * lon (lon) float32 0.0 5.0 10.0 15.0 20.0 ... 340.0 345.0 350.0 355.0 * lat (lat) float32 -90.0 -86.0 -82.0 -78.0 -74.0 ... 78.0 82.0 86.0 90.0 time object 0030-02-01 00:00:00 Attributes: long_name: area weights normalized by max units:
arwt_norm_max.plot(robust=True);
Needed some help getting scalar output from 2D weighted means function: https://github.com/NCAR/esmlab/issues/28
This is the purpose of .load().values
TSA_glob_mean_arwt_sum = esmlab.statistics.weighted_mean(TSA_amean, dim=['lat','lon'], weights=arwt_norm_sum).load().values
TSA_glob_mean_arwt_sum
array(283.07242, dtype=float32)
TSA_glob_mean_arwt_max = esmlab.statistics.weighted_mean(TSA_amean, dim=['lat','lon'], weights=arwt_norm_max).load().values
TSA_glob_mean_arwt_max
array(283.07245, dtype=float32)
TSA_glob_mean_arwt = esmlab.statistics.weighted_mean(TSA_amean, dim=['lat','lon'], weights=arwt).load().values
TSA_glob_mean_arwt
array(283.07245, dtype=float32)
Result is the same as normalizing by the maximum
TSA_glob_mean_arwt_nolf = esmlab.statistics.weighted_mean(TSA_amean, dim=['lat','lon'], weights=area).load().values
TSA_glob_mean_arwt_nolf
array(284.36517, dtype=float32)
lats = fils_in.lat
# np.cos expects input in radians
cwt = np.cos(lats * np.pi / 180.)
# NOTE: this doesn't account for landfrac...but how to do that when cwt is indexed only by latitude?
TSA_glob_mean_cwt = esmlab.statistics.weighted_mean(TSA_amean, dim=['lat'], weights=cwt).mean().values
TSA_glob_mean_cwt
array(279.76453, dtype=float32)
Most straightforward solution is probably to figure out how to write or convert the NCL gaus
function to python:
https://www.ncl.ucar.edu/Document/Functions/Built-in/gaus.shtml
This is based on the python tutorial xarray notebook:
https://github.com/ncar-hackathons/hands-on-examples/blob/master/scientific-computing/xarray.ipynb
The area element for lat-lon coordinates is
$$ \delta A = R^2 \delta \phi \delta \lambda \cos(\phi) $$where $\phi$ is latitude, $\delta \phi$ is the spacing of the points in latitude, $\delta \lambda$ is the spacing of the points in longitude, and $R$ is Earth's radius. (In this formula, $\phi$ and $\lambda$ are measured in radians.) Let's use xarray to create the weight factor.
R = 6.37e6
# we know already that the spacing of the points is 4 degrees latitude, 5 degrees longitude
dϕ = np.deg2rad(4.)
dλ = np.deg2rad(5.)
dA = R**2 * dϕ * dλ * np.cos(np.deg2rad(lats))
dA.plot();
dA.where(TSA_amean.notnull())
<xarray.DataArray 'lat' (lat: 46, lon: 72)> dask.array<where, shape=(46, 72), dtype=float32, chunksize=(46, 72), chunktype=numpy.ndarray> Coordinates: * lat (lat) float32 -90.0 -86.0 -82.0 -78.0 -74.0 ... 78.0 82.0 86.0 90.0 * lon (lon) float32 0.0 5.0 10.0 15.0 20.0 ... 340.0 345.0 350.0 355.0
pixel_area = dA.where(TSA_amean.notnull())
pixel_area.plot(robust=True);
# why are there negative values in the colorbar?
# NOTE: this also doesn't account for landfrac...
total_land_area = pixel_area.sum(dim=['lat', 'lon'])
TSA_glob_mean_ae = ((TSA_amean * pixel_area).sum(dim=['lat', 'lon']) / total_land_area).values
TSA_glob_mean_ae
array(284.4086, dtype=float32)
Try again, accounting for landfrac
pixel_area_lf = dA.where(TSA_amean.notnull())*landfrac
pixel_area_lf.plot(robust=True);
total_land_area_lf = pixel_area_lf.sum(dim=['lat', 'lon'])
TSA_glob_mean_ae_lf = ((TSA_amean * pixel_area_lf).sum(dim=['lat', 'lon']) / total_land_area_lf).values
TSA_glob_mean_ae_lf
array(283.12924, dtype=float32)
Now the result is closer to arwt methods
This function is equivalent to esmlab weighted mean, and result will vary similarly depending on choice of weights
def horizontal_weighted_mean(var, wgts):
return (var * wgts).sum(dim=['lat', 'lon']) / wgts.sum(dim=['lat', 'lon'])
TSA_glob_mean_geocat_sum = horizontal_weighted_mean(TSA_amean, arwt_norm_sum).values
TSA_glob_mean_geocat_sum
array(283.07242, dtype=float32)
TSA_glob_mean_geocat_max = horizontal_weighted_mean(TSA_amean, arwt_norm_max).values
TSA_glob_mean_geocat_max
array(283.07245, dtype=float32)
Without normalizing weights
TSA_glob_mean_geocat = horizontal_weighted_mean(TSA_amean, arwt).values
TSA_glob_mean_geocat
array(283.07245, dtype=float32)
Without landfrac
TSA_glob_mean_geocat_nolf = horizontal_weighted_mean(TSA_amean, area).values
TSA_glob_mean_geocat_nolf
array(284.36517, dtype=float32)
fig, ax = plt.subplots(figsize=(12,6))
means = [TSA_glob_mean_arwt_sum, TSA_glob_mean_arwt_max, TSA_glob_mean_arwt, TSA_glob_mean_arwt_nolf, TSA_glob_mean_cwt, TSA_glob_mean_ae, TSA_glob_mean_ae_lf, TSA_glob_mean_geocat]
x = np.arange(len(means))
means_c = [m - 273.15 for m in means]
# color code bars based on value - though hard to distinguish close values
from matplotlib import cm
color_vals = np.array(means_c)
colors = cm.hsv(color_vals / float(max(color_vals)))
bar_plot = plt.bar(x,means_c,color = colors)
plt.title("Comparing Global Mean Methods")
plt.ylabel("2m air temperature ($\degree$C)")
plt.xticks(x, ('arwt sum','arwt max', 'arwt no norm', 'area no lf', 'cos lat', 'area element', 'area element lf', 'geocat'))
# add text labels
labels_plot = [round(m, 5) for m in means_c] # round to 4 decimal places
def autolabel(rects):
for idx,rect in enumerate(bar_plot):
height = rect.get_height()
ax.text(rect.get_x() + rect.get_width()/2., 1.02*height,
labels_plot[idx],
ha='center', va='bottom', rotation=0)
autolabel(bar_plot)
plt.ylim(0,13);
Sidenote: global mean TSA values may seem low, but these are land only averages
# Another way to plot: points instead of bars
plt.figure(figsize=(12,6))
plt.plot(means_c, 'mo')
plt.title("Comparing Global Mean Methods")
plt.ylabel("2m air temperature ($\degree$C)")
plt.xticks(x, ('arwt sum','arwt max','arwt no norm','area no lf', 'cos lat', 'area element', 'area element lf', 'geocat'));
# Most of CONUS
lat_min = 30
lat_max = 50
lon_min = 230
lon_max = 290
TSA_amean_US = TSA_amean.sel(lat=slice(lat_min,lat_max), lon=slice(lon_min,lon_max))
TSA_amean_US
<xarray.DataArray 'TSA' (lat: 6, lon: 13)> dask.array<getitem, shape=(6, 13), dtype=float32, chunksize=(6, 13), chunktype=numpy.ndarray> Coordinates: * lon (lon) float32 230.0 235.0 240.0 245.0 ... 275.0 280.0 285.0 290.0 * lat (lat) float32 30.0 34.0 38.0 42.0 46.0 50.0 Attributes: long_name: 2m air temperature units: K cell_methods: time: mean
TSA_US_mean_arwt_sum = esmlab.statistics.weighted_mean(TSA_amean_US, dim=['lat','lon'], weights=arwt_norm_sum.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).load().values
TSA_US_mean_arwt_sum
array(283.41443, dtype=float32)
TSA_US_mean_arwt_max = esmlab.statistics.weighted_mean(TSA_amean_US, dim=['lat','lon'], weights=arwt_norm_max.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).load().values
TSA_US_mean_arwt_max
array(283.4144, dtype=float32)
TSA_US_mean_arwt = esmlab.statistics.weighted_mean(TSA_amean_US, dim=['lat','lon'], weights=arwt.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).load().values
TSA_US_mean_arwt
array(283.41446, dtype=float32)
TSA_US_mean_arwt_nolf = esmlab.statistics.weighted_mean(TSA_amean_US, dim=['lat','lon'], weights=area.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).load().values
TSA_US_mean_arwt_nolf
array(284.4153, dtype=float32)
TSA_US_mean_cwt = esmlab.statistics.weighted_mean(TSA_amean_US, dim=['lat'], weights=cwt.sel(lat=slice(lat_min,lat_max))).mean().values
TSA_US_mean_cwt
array(284.19653, dtype=float32)
CONUS_pixel_area = pixel_area.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))
CONUS_pixel_area.plot(robust=True)
<matplotlib.collections.QuadMesh at 0x7f4926100fd0>
CONUS_land_area = CONUS_pixel_area.sum(dim=['lat', 'lon'])
TSA_US_mean_ae = ((TSA_amean_US * CONUS_pixel_area).sum(dim=['lat', 'lon']) / CONUS_land_area).values
TSA_US_mean_ae
array(284.4153, dtype=float32)
CONUS_pixel_area_lf = CONUS_pixel_area*landfrac.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))
CONUS_pixel_area_lf.plot(robust=True);
CONUS_land_area_lf = CONUS_pixel_area_lf.sum(dim=['lat', 'lon'])
TSA_US_mean_ae_lf = ((TSA_amean_US * CONUS_pixel_area_lf).sum(dim=['lat', 'lon']) / CONUS_land_area_lf).values
TSA_US_mean_ae_lf
array(283.4145, dtype=float32)
TSA_US_mean_geocat_sum = horizontal_weighted_mean(TSA_amean_US, arwt_norm_sum.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).values
TSA_US_mean_geocat_sum
array(283.41443, dtype=float32)
TSA_US_mean_geocat_max = horizontal_weighted_mean(TSA_amean_US, arwt_norm_max.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).values
TSA_US_mean_geocat_max
array(283.4144, dtype=float32)
TSA_US_mean_geocat = horizontal_weighted_mean(TSA_amean_US, arwt.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).values
TSA_US_mean_geocat
array(283.41446, dtype=float32)
TSA_US_mean_geocat_nolf = horizontal_weighted_mean(TSA_amean_US, area.sel(lat=slice(lat_min,lat_max),lon=slice(lon_min,lon_max))).values
TSA_US_mean_geocat_nolf
array(284.4153, dtype=float32)
fig, ax = plt.subplots(figsize=(12,6))
means = [TSA_US_mean_arwt_sum, TSA_US_mean_arwt_max, TSA_US_mean_arwt, TSA_US_mean_arwt_nolf, TSA_US_mean_cwt, TSA_US_mean_ae, TSA_US_mean_ae_lf, TSA_US_mean_geocat]
x = np.arange(len(means))
means_c = [m - 273.15 for m in means]
# color code bars based on value - though hard to distinguish close values
from matplotlib import cm
color_vals = np.array(means_c)
colors = cm.hsv(color_vals / float(max(color_vals)))
bar_plot = plt.bar(x,means_c,color = colors)
plt.title("Comparing Regional Mean Methods")
plt.ylabel("CONUS 2m air temperature ($\degree$C)")
plt.xticks(x, ('arwt sum','arwt max', 'arwt no norm', 'area no lf', 'cos lat', 'area element', 'area element lf', 'geocat'))
# add text labels
labels_plot = [round(m, 5) for m in means_c] # round to 4 decimal places
def autolabel(rects):
for idx,rect in enumerate(bar_plot):
height = rect.get_height()
ax.text(rect.get_x() + rect.get_width()/2., 1.02*height,
labels_plot[idx],
ha='center', va='bottom', rotation=0)
autolabel(bar_plot)
plt.ylim(0,13);
Interesting, this time cos lat biases high, versus globally it biases low
Everything else is similar to global relationships