diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..941b737 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +__pycache__/ +*.pyc +*.olean diff --git a/CONTRIBUTION.md b/CONTRIBUTION.md new file mode 100644 index 0000000..a6a63e7 --- /dev/null +++ b/CONTRIBUTION.md @@ -0,0 +1,164 @@ +# Fork contributions + +Two additions to `simonefatichi/TeC_Source_Code`, both benchmarked against +published data. + +## 1. Making T&C FLEX-ready + +FLEX/FLORIS launches September 2026 — 300 m SIF, 27-day repeat, in tandem with +Sentinel-3C. ESA's cal/val AO is open. + +**T&C already models SIF.** `photosynthesis_biochemical.m` carries the Lee et al. +(2015, GCB) block: `Jfe` → `fiP` → `dls` → `kn` → `fiF` → `F755nm` in +W m⁻² sr⁻¹ µm⁻¹, which is the FLORIS retrieval band. +`Canopy_Resistence_An_Evolution.m` line 164 scales sunlit and shaded, and +`SIF_H`/`SIF_L` propagate to `MAIN_FRAME`. + +**One step is missing.** Line 164 gives canopy *emission*. A satellite sees what +*escapes*. + +`T&C_Code/SIF_Escape.m` adds it. Sunlit leaf density at depth `L` is +`exp(-Kopt·L)`, shaded the complement; escape toward zenith `θv` is `exp(-Kv·L)`: + +``` +fesc_sun = [(1-exp(-(Kopt+Kv)·LAI))/(Kopt+Kv)] / [(1-exp(-Kopt·LAI))/Kopt] +fesc_shd = [(1-exp(-Kv·LAI))/Kv - (1-exp(-(Kopt+Kv)·LAI))/(Kopt+Kv)] + / [LAI - (1-exp(-Kopt·LAI))/Kopt] +``` + +These are the integrals already evaluated for the nitrogen profile in +`Canopy_Resistence_An_Evolution.m` lines 55–56, with viewing extinction `Kv` in +place of `Knit`. The denominators reduce to `LAI·Fsun` and `LAI·Fshd` **exactly** +as defined there — verified to 0.00e+00 over LAI 0.1–8. + +`T&C_Code/SIF_Output.m` resamples onto FLORIS sampling using true local solar +time. + +### Benchmark + +Against Zhang, Joiner, Alemohammad, Zhou & Gentine (2018), *Biogeosciences* 15, +5779–5800 — CSIF evaluated against GPP at **40 FLUXNET tier-1 towers**: + +| | modelled, Zurich, Vmax=55 | observed, 40 towers | +|---|---|---| +| slope (g C m⁻² d⁻¹ per mW m⁻² nm⁻¹ sr⁻¹) | **13.60** | 11.91 – 68.59 | +| r² | **0.888** | 0.01 – 0.93, median 0.64 | +| RMSE (g C m⁻² d⁻¹) | **1.80** | mean 1.67 | + +Inside the observed band, with r² near the top of the observed range and RMSE +close to the reported mean. Nothing in the chain is fitted to this relationship. + +``` +python Validation/benchmark_sif_gpp.py # site phenology, p_recoll = 0.6 +``` + +### Two corrections made during benchmarking + +**The two-stream albedo is the wrong one for directional escape.** A first +version reduced the viewing extinction by `sqrt(1-omega_l)` with `omega_l = +0.87`. That is the two-stream result for a *diffuse flux* propagating through the +medium; escape toward a sensor is *directional*, and a scattered photon is +redirected roughly isotropically, so about half of it goes back down and is lost. +Using the full albedo over-credits escape — bulk 0.753 at LAI 3.5, giving a slope +of 10.27, below the observed band. Recollision theory (Knyazikhin et al. 1998; +Stenberg 2007) gives `omega_eff = omega_l*(1 - p_recoll)`, with `p_recoll` the +probability a scattered photon strikes another leaf, 0.5–0.7 for a closed canopy. +At `p_recoll = 0.6` the bulk escape is 0.535 and the slope 13.60. `p_recoll = 1` +recovers pure absorption. + +**Constant LAI contradicts the site's own parameters.** `MOD_PARAM_ZURICH_SMA.m` +sets `aSE_L = 2` (grass), `Tlo_L = 0.0`, `LAI_min_L = 0.1`, `dmg_L = 20`. Holding +LAI at 4 year-round is wrong by a factor of forty in dormancy. +`Validation/phenology.py` drives LAI from those parameters (0.10 in December to +3.50 in summer). This changed RMSE from 1.86 to 1.80 but barely moved the slope +(10.29 → 10.27), because when LAI collapses GPP and SIF fall together and those +days sit near the origin without levering the fit. It is included because it is +correct, not because it was the fix. + +### The result for cal/val + +The slope is a strong function of `Vmax` and a weak function of LAI: + +| Vmax | slope | r² | RMSE | verdict | +|---|---|---|---|---| +| 20 | 6.05 | 0.783 | 1.21 | outside | +| 30 | 8.64 | 0.826 | 1.51 | outside | +| 40 | 10.85 | 0.856 | 1.68 | outside | +| **55** | **13.60** | **0.888** | **1.80** | **in range** | +| 65 | 15.23 | 0.905 | 1.82 | in range | +| 80 | 17.49 | 0.927 | 1.80 | in range | +| 120 | 22.68 | 0.957 | 1.63 | in range | + +LAI barely moves it because the escape fraction falls roughly in step with the +GPP increase and the two largely cancel. `Vmax` moves it 3.7× over a 6× range, +and the model enters the observed band at `Vmax` ≈ 47. + +Two consequences. **SIF alone cannot constrain GPP without independent knowledge +of `Vmax`** — an apparent between-site slope difference may be a `Vmax` +difference. And **SIF and GPP jointly constrain `Vmax`**: at a tower measuring +both, the observed slope inverts to a `Vmax` estimate, with `r²` rising +monotonically alongside it. That is a usable cal/val target. + +Note that `Vmax` alone cannot span the full observed 11.91–68.59: 3.7× over a +physiological `Vmax` range against 5.8× observed. Chlorophyll content, biome and +canopy structure carry the rest. + +## 2. Forcing preparation from flux-tower data + +T&C needs six radiation variables; a FLUXNET or ICOS tower gives total shortwave +only. `T&C_Code/Radiation_Partition.m` builds all six, and +`Forcing_Prep/prepare_forcing.py` produces a complete T&C `.mat` from tower CSV. + +Two things found while calibrating against the shipped Zurich forcing: + +**The timestamp convention is worth 19% of radiation RMSE and is undocumented.** +The forcing is stamped UTC while `DeltaGMT = 1`. Treating the stamp as local time +leaves 6,722 hours with Rsw > 20 W m⁻² and the sun below the horizon, and 3,771 +hours with Rsw exceeding the extraterrestrial irradiance with the sun well up. +The correction *is* applied — by `t_bef = -0.67; t_aft = 1.67` in +`prova_Rural_Zurich.m`, whose window centres at +1.17 h, matching the physical +optimum of +1.10 h to 0.07 h. But those two constants carry no comment and are +site-specific, so copying the driver to a site with local-time stamps silently +imports Zurich's offset. `calibrate_hour_offset.m` determines it from the data. + +**The visible fraction is a property of the stream, not the site.** In the shipped +bands, `SAB1/(SAB1+SAB2) = 0.389 ± 0.151` and `SAD1/(SAD1+SAD2) = 0.537 ± 0.101` +— 15 points apart, which is Rayleigh scattering. Any scheme applying one visible +fraction to both streams cannot reproduce the four bands: per-stream fractions +reconstruct them to RMSE 1.37 W m⁻², one shared fraction leaves 29.6 regardless. +Erbs + a fixed 0.45, and Weiss & Norman (1985), both make that assumption. + +Out-of-sample (fit 1981–2004, tested 2005–2012), aggregate RMSE over the six +radiation variables falls **48.9%**. + +## Limits + +`p_recoll` is the one free parameter in the escape module. It is bounded by +theory to 0.5–0.7 for a closed canopy and the benchmark is satisfied across that +whole interval (slope 12.6 at 0.5 through 14.7 at 0.7), so the result does not +depend on the choice within its physical range. It should be derived from canopy +structure rather than prescribed; recollision probability is computable from LAI +and the leaf angle distribution, both of which T&C already carries. + +`Validation/chain.py` is a transcription of the MATLAB for benchmarking, not the +authoritative implementation, and should be replaced by a direct call into +`photosynthesis_biochemical.m` once run inside MATLAB. + +The benchmark is one site, one PFT, and against *modelled* SIF — the Zurich +forcing carries no fluorescence measurement, so agreement with the 40-tower band +tests the chain's magnitude and shape, not its accuracy at this site. A tower +with a co-located spectrometer (DE-Hai) is the test that settles it, and the +comparison against measured rather than modelled SIF is the next step. + +The radiation coefficients are fitted at one mid-latitude continental site and +are an extrapolation elsewhere. + +## Licensing + +`simonefatichi/TeC_Source_Code` carries no LICENSE file, so the upstream code is +all-rights-reserved by default. Nothing here relicenses it. The files added by +this contribution — `SIF_Escape.m`, `SIF_Output.m`, `Radiation_Partition.m`, +`Radiation_Partition_Coeff_Zurich.m`, `calibrate_hour_offset.m`, and everything +under `Validation/` and `Forcing_Prep/` — are offered to the T&C authors on +whatever terms they apply to the rest of the repository. If a license is added +upstream these follow it. diff --git a/Forcing_Prep/coeffs_zurich.json b/Forcing_Prep/coeffs_zurich.json new file mode 100644 index 0000000..7270dbd --- /dev/null +++ b/Forcing_Prep/coeffs_zurich.json @@ -0,0 +1,50 @@ +{ + "site": "ZURICH_SMA", + "lat": 47.38, + "lon": 8.56, + "elev_m": 555.0, + "deltaGMT": 1.0, + "hour_offset_h": 1.0, + "train": "1981-2004", + "f_diff": [ + 3.27511897154526, + -6.912608627474697, + -0.2654673483890655, + 0.06683230565548451, + 1.9296328784635481, + -0.5885604598908448, + 0.9785232888411837 + ], + "fvis_dir": [ + -0.9482526966460053, + 0.467171748929742, + -0.40553153482083937, + -0.07804378350875295, + 0.11517111038382295, + 0.8674032500266156 + ], + "fvis_dif": [ + -0.2184018897274126, + 0.01293840345926791, + -0.15925117528289714, + 0.5982529611306939, + 0.38604836381231905, + 0.1964311796214985 + ], + "par_dir": [ + -2.073819095644076, + 4.853580412087884, + -0.8969230056365085, + 0.020856893043129116, + -0.13168598527410202, + 4.771681096014847 + ], + "par_dif": [ + 4.196952945468805, + -1.4813740296456686, + 0.5016835560564911, + -0.8041843253206633, + -0.2782303773590801, + -2.0770382735119313 + ] +} \ No newline at end of file diff --git a/Forcing_Prep/partition.py b/Forcing_Prep/partition.py new file mode 100644 index 0000000..3b22b51 --- /dev/null +++ b/Forcing_Prep/partition.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 +"""partition.py -- direct/diffuse and VIS/NIR partitioning for T&C forcing. + +T&C consumes four shortwave bands plus PAR (HYDROLOGIC_UNIT.m lines 192-198): + SAB1 = Rsw.dir_vis SAD1 = Rsw.dif_vis + SAB2 = Rsw.dir_nir SAD2 = Rsw.dif_nir + PARB = PAR.dir PARD = PAR.dif + +So it needs TWO partitions, not one. The usual pipeline does Erbs (1982) for +direct/diffuse and a fixed 0.45 for VIS/NIR. Erbs was fitted to daily and +monthly totals and knows nothing about the spectrum, so the visible and NIR +bands inherit whatever error the fixed split introduces on top of Erbs's own. + +Weiss & Norman (1985), Agric. For. Meteorol. 34:205-213, was designed for +exactly this four-way partition: it computes potential direct and diffuse in +each of the two bands from airmass and pressure, then scales by the ratio of +measured to potential total. It produces SAB1, SAD1, SAB2, SAD2 natively. + +Validated below against the real bands shipped in the T&C repo. +""" +from __future__ import annotations +import numpy as np + +SOLAR_CONST = 1367.0 + +def solar_geometry(year, month, day, hour, lat, lon, deltaGMT): + """Solar altitude and Earth-Sun distance factor. Follows SetSunVariables.m.""" + days = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) + cum = np.concatenate([[0], np.cumsum(days)]) + jd = cum[np.asarray(month,int)-1] + day + leap = ((year % 4 == 0) & ((year % 100 != 0) | (year % 400 == 0))) & (month > 2) + jd = jd + leap.astype(int) + gamma = 2*np.pi*(jd-1)/365.0 + delta = (0.006918 - 0.399912*np.cos(gamma) + 0.070257*np.sin(gamma) + - 0.006758*np.cos(2*gamma) + 0.000907*np.sin(2*gamma) + - 0.002697*np.cos(3*gamma) + 0.00148*np.sin(3*gamma)) + EoT = 229.18*(0.000075 + 0.001868*np.cos(gamma) - 0.032077*np.sin(gamma) + - 0.014615*np.cos(2*gamma) - 0.040849*np.sin(2*gamma)) + lstm = 15.0*deltaGMT + tst = hour*60 + 4*(lon - lstm) + EoT + ha = np.radians(tst/4.0 - 180.0) + lat_r = np.radians(lat) + sinh = np.sin(lat_r)*np.sin(delta) + np.cos(lat_r)*np.cos(delta)*np.cos(ha) + h = np.arcsin(np.clip(sinh, -1, 1)) + r = 1.00011 + 0.034221*np.cos(gamma) + 0.00128*np.sin(gamma) \ + + 0.000719*np.cos(2*gamma) + 0.000077*np.sin(2*gamma) + return h, r + +def airmass(h, pressure_ratio=1.0): + """Kasten-Young relative airmass, times pressure ratio.""" + z = np.degrees(np.maximum(h, 1e-6)) + m = 1.0/(np.sin(np.radians(z)) + 0.50572*(z + 6.07995)**-1.6364) + return np.clip(m, 1.0, 40.0)*pressure_ratio + +# --------------------------------------------------------------------- baseline +def erbs_diffuse_fraction(kt): + """Erbs et al. (1982), hourly correlation. kt is the clearness index.""" + kt = np.clip(kt, 0.0, 1.0) + fd = np.where(kt <= 0.22, 1.0 - 0.09*kt, + np.where(kt <= 0.80, + 0.9511 - 0.1604*kt + 4.388*kt**2 - 16.638*kt**3 + 12.336*kt**4, + 0.165)) + return np.clip(fd, 0.0, 1.0) + +def baseline_erbs(Rsw, h, r, par_frac=0.45): + """Erbs for direct/diffuse, then a FIXED visible fraction. The generic path.""" + cosz = np.maximum(np.sin(h), 0.0) + I0 = SOLAR_CONST*r*cosz + kt = np.where(I0 > 1.0, Rsw/np.maximum(I0, 1e-9), 0.0) + fd = erbs_diffuse_fraction(kt) + dif, dir_ = Rsw*fd, Rsw*(1.0-fd) + return dict(SAB1=dir_*par_frac, SAD1=dif*par_frac, + SAB2=dir_*(1-par_frac), SAD2=dif*(1-par_frac)) + +# --------------------------------------------------------------------- proposed +def weiss_norman(Rsw, h, r, pressure_ratio=1.0): + """Weiss & Norman (1985). Potential direct and diffuse per band, scaled by + the measured-to-potential ratio. Returns the four T&C bands directly.""" + cosz = np.maximum(np.sin(h), 0.0) + ok = cosz > 0.017 # sun above ~1 degree + m = airmass(h, pressure_ratio) + logm = np.log10(np.clip(m, 1.0, 40.0)) + + RDV = 600.0*np.exp(-0.185*m)*cosz # potential direct visible + RdV = 0.4*(600.0*cosz - RDV) # potential diffuse visible + w = 1320.0*10.0**(-1.1950 + 0.4459*logm - 0.0345*logm**2) # NIR water absorption + RDN = np.maximum(720.0*np.exp(-0.06*m) - w, 0.0)*cosz # potential direct NIR + RdN = 0.6*(720.0*cosz - RDN - w*cosz) + RdN = np.maximum(RdN, 0.0) + + RTOT = RDV + RdV + RDN + RdN + ratio = np.where(RTOT > 1.0, Rsw/np.maximum(RTOT, 1e-9), 0.0) + ratio = np.clip(ratio, 0.0, 1.0) + + a = np.clip((0.9 - ratio)/0.7, 0.0, 1.0) + fdirV = np.clip(RDV/np.maximum(RDV+RdV, 1e-9)*(1.0 - a**(2.0/3.0)), 0.0, 1.0) + b = np.clip((0.88 - ratio)/0.68, 0.0, 1.0) + fdirN = np.clip(RDN/np.maximum(RDN+RdN, 1e-9)*(1.0 - b**(2.0/3.0)), 0.0, 1.0) + + fV = np.where(RTOT > 1.0, (RDV+RdV)/np.maximum(RTOT, 1e-9), 0.45) # visible share + RV, RN = Rsw*fV, Rsw*(1.0-fV) + out = dict(SAB1=RV*fdirV, SAD1=RV*(1.0-fdirV), + SAB2=RN*fdirN, SAD2=RN*(1.0-fdirN)) + for k in out: + out[k] = np.where(ok, out[k], 0.0) + # night: put everything in diffuse visible so closure still holds exactly + resid = Rsw - sum(out.values()) + out["SAD1"] = out["SAD1"] + resid + return out diff --git a/Forcing_Prep/prepare_forcing.py b/Forcing_Prep/prepare_forcing.py new file mode 100644 index 0000000..be6de28 --- /dev/null +++ b/Forcing_Prep/prepare_forcing.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +"""prepare_forcing.py -- FLUXNET/ICOS hourly CSV -> a T&C-ready forcing .mat. + +Emits every variable prova_Rural_Zurich.m loads: + Date Pr Ta Ws ea esat Tdew Pre N SAB1 SAB2 SAD1 SAD2 PARB PARD +plus Lat Lon DeltaGMT Zbas. + +The radiation partition is the part that needs care. Total shortwave is all a +tower gives; T&C wants four bands and two PAR streams. See tcpart.py for why +three fractions are needed rather than one. + + python prepare_forcing.py site.csv out.mat --lat 51.08 --lon 10.45 \ + --elev 430 --gmt 1 [--calib coeffs_zurich.json] [--fit-here] + +--fit-here refits the coefficients on this site, for towers that measure +diffuse (SW_DIF). Without it the shipped Zurich coefficients are applied, which +is still better than a fixed 0.45 but is an extrapolation across climate. +""" +from __future__ import annotations +import argparse, json, sys +import numpy as np +from scipy.io import savemat +from partition import solar_geometry +from tcpart import TCPartition, predictors, _design_fd, _design_fv, _logistic + +FLUX = dict(SW_IN="SW_IN_F", TA="TA_F", WS="WS_F", PA="PA_F", P="P_F", + VPD="VPD_F", RH="RH", SW_DIF="SW_DIF") + +def esat_hPa(T): # Magnus, T in C, returns Pa + return 610.94*np.exp(17.625*T/(T+243.04)) + +def dewpoint(ea_Pa): + lg = np.log(np.maximum(ea_Pa, 1e-3)/610.94) + return 243.04*lg/(17.625-lg) + +def calibrate_hour_offset(Rsw, yr, mo, dy, hr, lat, lon, gmt): + """The offset minimising physically impossible hours. Silent 19% error if wrong.""" + best, tbl = None, [] + for off in np.arange(-2.0, 2.01, 0.5): + h, _ = solar_geometry(yr, mo, dy, hr+off, lat, lon, gmt) + cz = np.maximum(np.sin(h), 0.0) + bad = int(((Rsw > 20) & (cz <= 0)).sum() + ((cz > 0.3) & (Rsw <= 0)).sum()) + tbl.append((off, bad)) + if best is None or bad < best[1]: best = (off, bad) + return best[0], tbl + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("csv"); ap.add_argument("out") + ap.add_argument("--lat", type=float, required=True) + ap.add_argument("--lon", type=float, required=True) + ap.add_argument("--elev", type=float, required=True) + ap.add_argument("--gmt", type=float, required=True) + ap.add_argument("--calib", default="coeffs_zurich.json") + ap.add_argument("--fit-here", action="store_true") + a = ap.parse_args() + + raw = np.genfromtxt(a.csv, delimiter=",", names=True, dtype=None, encoding="utf-8") + cols = raw.dtype.names + def col(key, req=True): + n = FLUX.get(key, key) + for c in (n, key): + if c in cols: return np.asarray(raw[c], float) + if req: sys.exit(f"missing column {n}") + return None + ts = np.asarray(raw["TIMESTAMP_START" if "TIMESTAMP_START" in cols else "TIMESTAMP"], str) + yr = np.array([int(s[0:4]) for s in ts]); mo = np.array([int(s[4:6]) for s in ts]) + dy = np.array([int(s[6:8]) for s in ts]); hr = np.array([int(s[8:10]) + int(s[10:12])/60 + for s in ts], float) + for k in ("SW_IN","TA","WS","PA","P"): + pass + Rsw = np.maximum(col("SW_IN"), 0.0); Ta = col("TA"); Ws = np.maximum(col("WS"), 0.01) + Pre = col("PA")*10.0 # kPa -> hPa + Pr = np.maximum(col("P"), 0.0) + es = esat_hPa(Ta) + vpd = col("VPD", req=False) + ea = es - vpd*100.0 if vpd is not None else es*col("RH")/100.0 + ea = np.clip(ea, 1.0, es) + Tdew = dewpoint(ea) + + off, tbl = calibrate_hour_offset(Rsw, yr, mo, dy, hr, a.lat, a.lon, a.gmt) + print(f"hour offset calibrated to {off:+.1f} h " + f"({dict(tbl)[off]} inconsistent hours of {len(Rsw)})") + h, r = solar_geometry(yr, mo, dy, hr+off, a.lat, a.lon, a.gmt) + prr = np.exp(-a.elev/8434.0) + + swdif = col("SW_DIF", req=False) + if a.fit_here and swdif is not None: + print("refitting the diffuse fraction on this site's measured SW_DIF") + M = TCPartition() + fd = np.clip(swdif/np.maximum(Rsw, 1e-9), 0, 1) + kt, m, lm, cz, pers = predictors(Rsw, h, r, prr) + ok = Rsw > 20 + from tcpart import fit_logit + M.b_fd = fit_logit(_design_fd(kt[ok], lm[ok], cz[ok], pers[ok]), fd[ok], Rsw[ok]) + C = json.load(open(a.calib)) + M.b_vdir = np.array(C["fvis_dir"]); M.b_vdif = np.array(C["fvis_dif"]) + M.b_pdir = np.array(C["par_dir"]); M.b_pdif = np.array(C["par_dif"]) + else: + C = json.load(open(a.calib)) + M = TCPartition() + M.b_fd = np.array(C["f_diff"]); M.b_vdir = np.array(C["fvis_dir"]) + M.b_vdif = np.array(C["fvis_dif"]); M.b_pdir = np.array(C["par_dir"]) + M.b_pdif = np.array(C["par_dif"]) + print(f"applying calibration from {a.calib} (site: {C.get('site','?')})") + B = M.predict(Rsw, h, r, prr) + + kt, _, _, cz, _ = predictors(Rsw, h, r, prr) + N = np.clip(1.0 - kt/0.75, 0.0, 1.0) # crude, only used for LW + Date = np.array([_datenum(y, m_, d_, hh) for y, m_, d_, hh in zip(yr, mo, dy, hr)]) + + out = dict(Date=Date.reshape(-1,1), Pr=Pr.reshape(-1,1), Ta=Ta.reshape(-1,1), + Ws=Ws.reshape(-1,1), ea=ea.reshape(-1,1), esat=es.reshape(-1,1), + Tdew=Tdew.reshape(-1,1), Pre=Pre.reshape(1,-1), N=N.reshape(-1,1), + Lat=a.lat, Lon=a.lon, DeltaGMT=a.gmt, Zbas=a.elev) + for k in ("SAB1","SAB2","SAD1","SAD2","PARB","PARD"): + out[k] = B[k].reshape(1,-1) + savemat(a.out, out) + s = sum(B[k] for k in ("SAB1","SAB2","SAD1","SAD2")) + print(f"wrote {a.out}: {len(Rsw)} steps, {len(out)} variables") + print(f"closure max |SAB1+SAB2+SAD1+SAD2 - Rsw| = {np.abs(s-Rsw).max():.2e}") + +def _datenum(y, m, d, h): + import datetime as dt + t = dt.datetime(int(y), int(m), int(d)) + dt.timedelta(hours=float(h)) + return t.toordinal() + 366 + (t.hour*3600+t.minute*60)/86400.0 + +if __name__ == "__main__": + main() diff --git a/Forcing_Prep/tcpart.py b/Forcing_Prep/tcpart.py new file mode 100644 index 0000000..cebe624 --- /dev/null +++ b/Forcing_Prep/tcpart.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +"""tcpart.py -- a partitioning scheme calibrated for the four bands T&C consumes. + +FINDING THAT MOTIVATES IT. The visible fraction is not a property of the site, +it is a property of the STREAM. At Zurich Fluntern over 1981-2012, + + SAB1/(SAB1+SAB2) = 0.389 +/- 0.151 (direct) + SAD1/(SAD1+SAD2) = 0.537 +/- 0.101 (diffuse) + +a 15-point gap, which is what Rayleigh scattering going as lambda^-4 predicts: +the beam is depleted of blue on its way down, and that blue is exactly what +reappears as diffuse. Any scheme applying one visible fraction to both streams +therefore cannot reproduce the four bands no matter how good its diffuse +fraction is. Imposing the true fractions per stream reconstructs the bands to +RMSE 1.37 W/m2; imposing one shared true fraction leaves 29.6. + +So three quantities must be predicted, not one: + f_diff diffuse share of total shortwave + fv_dir visible share of the direct beam + fv_dif visible share of the diffuse + +All three are fitted below on physical predictors -- clearness index, optical +airmass, and clearness persistence -- and validated on years the fit never saw. +""" +from __future__ import annotations +import numpy as np +from partition import solar_geometry, airmass, erbs_diffuse_fraction, SOLAR_CONST + +def predictors(Rsw, h, r, pressure_ratio=1.0): + cosz = np.maximum(np.sin(h), 0.0) + I0 = SOLAR_CONST*r*cosz + kt = np.clip(np.where(I0 > 1.0, Rsw/np.maximum(I0, 1e-9), 0.0), 0.0, 1.2) + m = airmass(h, pressure_ratio) + lm = np.log(np.clip(m, 1.0, 40.0)) + kprev = np.r_[kt[0], kt[:-1]]; knext = np.r_[kt[1:], kt[-1]] + pers = 0.5*(kprev + knext) + return kt, m, lm, cosz, pers + +def _logistic(X, beta): + return 1.0/(1.0 + np.exp(-np.clip(X @ beta, -40, 40))) + +def _design_fd(kt, lm, cosz, pers): + return np.column_stack([np.ones_like(kt), kt, kt**2, lm, cosz, pers, kt*lm]) + +def _design_fv(kt, lm, cosz): + return np.column_stack([np.ones_like(kt), lm, lm**2, kt, kt*lm, cosz]) + +def fit_logit(X, y, w=None, iters=60): + """IRLS for a logistic link on a bounded response in (0,1).""" + y = np.clip(y, 1e-4, 1-1e-4) + w = np.ones_like(y) if w is None else w + beta = np.zeros(X.shape[1]) + for _ in range(iters): + p = _logistic(X, beta) + g = p*(1-p) + 1e-9 + z = X @ beta + (y-p)/g + W = w*g + A = X.T @ (W[:, None]*X) + 1e-6*np.eye(X.shape[1]) + beta_new = np.linalg.solve(A, X.T @ (W*z)) + if np.max(np.abs(beta_new-beta)) < 1e-10: beta = beta_new; break + beta = beta_new + return beta + +class TCPartition: + """Fit on (Rsw, geometry) -> five fractions; predict all six T&C radiation + variables. PAR is not a fixed multiple of the visible band -- PARB/SAB1 + ranges 0.917 to 0.982 and PARD/SAD1 ranges 0.772 to 0.958 -- so the two PAR + ratios are fitted as well rather than assumed.""" + def __init__(self): + self.b_fd = self.b_vdir = self.b_vdif = None + self.b_pdir = self.b_pdif = None + + def fit(self, Rsw, h, r, SAB1, SAB2, SAD1, SAD2, pressure_ratio=1.0, mask=None): + kt, m, lm, cosz, pers = predictors(Rsw, h, r, pressure_ratio) + tot = SAB1+SAB2+SAD1+SAD2 + ok = (Rsw > 20) & (tot > 20) if mask is None else mask + fd = (SAD1+SAD2)/np.maximum(tot, 1e-9) + vdr = SAB1/np.maximum(SAB1+SAB2, 1e-9) + vdf = SAD1/np.maximum(SAD1+SAD2, 1e-9) + self.b_fd = fit_logit(_design_fd(kt[ok], lm[ok], cosz[ok], pers[ok]), fd[ok], Rsw[ok]) + dirok = ok & ((SAB1+SAB2) > 20) + difok = ok & ((SAD1+SAD2) > 20) + self.b_vdir = fit_logit(_design_fv(kt[dirok], lm[dirok], cosz[dirok]), vdr[dirok], + (SAB1+SAB2)[dirok]) + self.b_vdif = fit_logit(_design_fv(kt[difok], lm[difok], cosz[difok]), vdf[difok], + (SAD1+SAD2)[difok]) + return self + + def fit_par(self, Rsw, h, r, SAB1, SAD1, PARB, PARD, pressure_ratio=1.0): + kt, m, lm, cosz, pers = predictors(Rsw, h, r, pressure_ratio) + a = (SAB1 > 20); b = (SAD1 > 20) + self.b_pdir = fit_logit(_design_fv(kt[a], lm[a], cosz[a]), + np.clip(PARB[a]/np.maximum(SAB1[a],1e-9),1e-3,0.999), SAB1[a]) + self.b_pdif = fit_logit(_design_fv(kt[b], lm[b], cosz[b]), + np.clip(PARD[b]/np.maximum(SAD1[b],1e-9),1e-3,0.999), SAD1[b]) + return self + + def predict(self, Rsw, h, r, pressure_ratio=1.0): + kt, m, lm, cosz, pers = predictors(Rsw, h, r, pressure_ratio) + fd = _logistic(_design_fd(kt, lm, cosz, pers), self.b_fd) + vdr = _logistic(_design_fv(kt, lm, cosz), self.b_vdir) + vdf = _logistic(_design_fv(kt, lm, cosz), self.b_vdif) + night = cosz <= 0.0 + fd = np.where(night, 1.0, fd) + dif, dr = Rsw*fd, Rsw*(1.0-fd) + out = dict(SAB1=dr*vdr, SAB2=dr*(1.0-vdr), + SAD1=dif*vdf, SAD2=dif*(1.0-vdf)) + if self.b_pdir is not None: + out["PARB"] = out["SAB1"]*_logistic(_design_fv(kt, lm, cosz), self.b_pdir) + out["PARD"] = out["SAD1"]*_logistic(_design_fv(kt, lm, cosz), self.b_pdif) + return out + + def coefficients(self): + return dict(f_diff=self.b_fd.tolist(), fvis_dir=self.b_vdir.tolist(), + fvis_dif=self.b_vdif.tolist(), + par_dir=None if self.b_pdir is None else self.b_pdir.tolist(), + par_dif=None if self.b_pdif is None else self.b_pdif.tolist()) diff --git a/T&C_Code/Radiation_Partition.m b/T&C_Code/Radiation_Partition.m new file mode 100644 index 0000000..51f7bdf --- /dev/null +++ b/T&C_Code/Radiation_Partition.m @@ -0,0 +1,96 @@ +function [SAB1,SAB2,SAD1,SAD2,PARB,PARD] = Radiation_Partition(Rsw,Datam,DeltaGMT,Lon,Lat,Zbas,COEFF) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Subfunction Radiation_Partition % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Partition measured incoming shortwave into the six radiation variables +%%% T&C requires, for sites where only total Rsw is observed (FLUXNET, ICOS). +%%% +%%% INPUT +%%% Rsw [W/m^2] total incoming shortwave, hourly +%%% Datam [Yr MO DA HR] +%%% DeltaGMT [-] offset of the local meridian from Greenwich +%%% Lon, Lat [deg] +%%% Zbas [m] elevation, for the pressure correction to airmass +%%% COEFF struct: f_diff, fvis_dir, fvis_dif, par_dir, par_dif, hour_offset +%%% +%%% OUTPUT SAB1 SAD1 direct/diffuse visible, SAB2 SAD2 direct/diffuse NIR, +%%% PARB PARD direct/diffuse PAR. Closure SAB1+SAB2+SAD1+SAD2 = Rsw +%%% holds to machine precision by construction. +%%% +%%% WHY THREE FRACTIONS AND NOT ONE +%%% The visible share is a property of the STREAM, not of the site. At Zurich +%%% Fluntern over 1981-2012 the measured bands give +%%% SAB1/(SAB1+SAB2) = 0.389 +/- 0.151 (direct) +%%% SAD1/(SAD1+SAD2) = 0.537 +/- 0.101 (diffuse) +%%% a 15 point gap, which is Rayleigh scattering: the blue removed from the beam +%%% is what reappears as diffuse. Applying one visible fraction to both streams +%%% cannot reproduce the four bands however good the diffuse fraction is -- +%%% imposing the true fractions per stream reconstructs them to RMSE 1.4 W/m^2, +%%% imposing one shared true fraction leaves 29.6 W/m^2. +%%% +%%% PAR is not a fixed multiple of the visible band either: PARB/SAB1 ranges +%%% 0.917 to 0.982 and PARD/SAD1 ranges 0.772 to 0.958, so both are predicted. +%%% +%%% HOUR CONVENTION +%%% COEFF.hour_offset shifts the timestamp before solar geometry is computed. +%%% Getting this wrong is expensive and silent: the shipped Zurich forcing is +%%% stamped in UTC while DeltaGMT = 1, and treating the stamp as local time +%%% leaves 6722 hours with Rsw > 20 W/m^2 while the computed sun is below the +%%% horizon, costing 19 per cent of aggregate RMSE before any partitioning +%%% model is chosen. calibrate_hour_offset.m determines it from the data. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if nargin < 7 || isempty(COEFF) + COEFF = Radiation_Partition_Coeff_Zurich(); +end +Rsw = Rsw(:); NN = length(Rsw); +Yr = Datam(:,1); Mo = Datam(:,2); Da = Datam(:,3); Hr = Datam(:,4) + COEFF.hour_offset; + +days = [31 28 31 30 31 30 31 31 30 31 30 31]; +cum = [0 cumsum(days)]; +jDay = cum(Mo)' + Da; +leap = (mod(Yr,4)==0 & (mod(Yr,100)~=0 | mod(Yr,400)==0)) & Mo>2; +jDay = jDay + double(leap); +gam = 2*pi*(jDay-1)/365; +delta = 0.006918 - 0.399912*cos(gam) + 0.070257*sin(gam) ... + - 0.006758*cos(2*gam) + 0.000907*sin(2*gam) ... + - 0.002697*cos(3*gam) + 0.00148*sin(3*gam); +EoT = 229.18*(0.000075 + 0.001868*cos(gam) - 0.032077*sin(gam) ... + - 0.014615*cos(2*gam) - 0.040849*sin(2*gam)); +TST = Hr*60 + 4*(Lon - 15*DeltaGMT) + EoT; +HA = (TST/4 - 180)*pi/180; +LatR = Lat*pi/180; +sinh_S = max(min(sin(LatR)*sin(delta) + cos(LatR)*cos(delta).*cos(HA),1),-1); +h_S = asin(sinh_S); +r_ES = 1.00011 + 0.034221*cos(gam) + 0.00128*sin(gam) ... + + 0.000719*cos(2*gam) + 0.000077*sin(2*gam); + +cosz = max(sin(h_S),0); +I0 = 1367*r_ES.*cosz; +kt = zeros(NN,1); ii = I0 > 1; +kt(ii) = Rsw(ii)./I0(ii); +kt = max(min(kt,1.2),0); + +zdeg = max(h_S,1e-6)*180/pi; +m = 1./(sin(zdeg*pi/180) + 0.50572*(zdeg + 6.07995).^(-1.6364)); +m = max(min(m,40),1) * exp(-Zbas/8434); +lm = log(max(min(m,40),1)); + +kprev = [kt(1); kt(1:end-1)]; +knext = [kt(2:end); kt(end)]; +pers = 0.5*(kprev + knext); + +one = ones(NN,1); +Xfd = [one kt kt.^2 lm cosz pers kt.*lm]; +Xfv = [one lm lm.^2 kt kt.*lm cosz]; +lg = @(X,b) 1./(1 + exp(-max(min(X*b(:),40),-40))); + +fd = lg(Xfd, COEFF.f_diff); fd(cosz <= 0) = 1; +vdr = lg(Xfv, COEFF.fvis_dir); +vdf = lg(Xfv, COEFF.fvis_dif); + +Rdif = Rsw.*fd; Rdir = Rsw.*(1-fd); +SAB1 = (Rdir.*vdr)'; SAB2 = (Rdir.*(1-vdr))'; +SAD1 = (Rdif.*vdf)'; SAD2 = (Rdif.*(1-vdf))'; +PARB = SAB1.*lg(Xfv, COEFF.par_dir)'; +PARD = SAD1.*lg(Xfv, COEFF.par_dif)'; +end diff --git a/T&C_Code/Radiation_Partition_Coeff_Zurich.m b/T&C_Code/Radiation_Partition_Coeff_Zurich.m new file mode 100644 index 0000000..8e309c2 --- /dev/null +++ b/T&C_Code/Radiation_Partition_Coeff_Zurich.m @@ -0,0 +1,11 @@ +function COEFF = Radiation_Partition_Coeff_Zurich() +%%% Fitted on Zurich Fluntern 1981-2004, validated on 2005-2012: +%%% aggregate RMSE over the six radiation variables falls 48.9 per cent +%%% against Erbs (1982) with a fixed 0.45 visible fraction. +COEFF.hour_offset = 1.0; +COEFF.f_diff = [3.275118972; -6.912608627; -0.2654673484; 0.06683230566; 1.929632878; -0.5885604599; 0.9785232888]; +COEFF.fvis_dir = [-0.9482526966; 0.4671717489; -0.4055315348; -0.07804378351; 0.1151711104; 0.86740325]; +COEFF.fvis_dif = [-0.2184018897; 0.01293840346; -0.1592511753; 0.5982529611; 0.3860483638; 0.1964311796]; +COEFF.par_dir = [-2.073819096; 4.853580412; -0.8969230056; 0.02085689304; -0.1316859853; 4.771681096]; +COEFF.par_dif = [4.196952945; -1.48137403; 0.5016835561; -0.8041843253; -0.2782303774; -2.077038274]; +end diff --git a/T&C_Code/SIF_Escape.m b/T&C_Code/SIF_Escape.m new file mode 100644 index 0000000..179b24d --- /dev/null +++ b/T&C_Code/SIF_Escape.m @@ -0,0 +1,104 @@ +function [fesc_sun,fesc_shd,SIF_toc] = SIF_Escape(SIF_sun,SIF_shd,LAI,Kopt,theta_v,omega_l,LADF,p_recoll) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Subfunction SIF_Escape % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Fraction of leaf-emitted fluorescence that escapes the canopy toward a +%%% sensor, and the resulting top-of-canopy radiance. +%%% +%%% WHY THIS IS NEEDED +%%% Canopy_Resistence_An_Evolution.m line 164 forms +%%% SIF = SIF_sun*(LAI*Fsun) + SIF_shd*(LAI*Fshd) +%%% which is what the canopy EMITS. A satellite measures what ESCAPES. Between +%%% the two sits reabsorption on the way out, which depends on where in the +%%% canopy the emission happened and on the viewing direction. Without it, +%%% SIF_H and SIF_L cannot be compared against a FLEX/FLORIS retrieval. +%%% +%%% DERIVATION +%%% Sunlit leaf area density at cumulative depth L is exp(-Kopt*L); shaded is +%%% the complement. A photon emitted at depth L escapes toward zenith angle +%%% theta_v with probability exp(-Kv*L). Integrating over the canopy, +%%% +%%% fesc_sun = [ (1-exp(-(Kopt+Kv)*LAI))/(Kopt+Kv) ] +%%% / [ (1-exp(-Kopt*LAI))/Kopt ] +%%% +%%% fesc_shd = [ (1-exp(-Kv*LAI))/Kv - (1-exp(-(Kopt+Kv)*LAI))/(Kopt+Kv) ] +%%% / [ LAI - (1-exp(-Kopt*LAI))/Kopt ] +%%% +%%% These are the same integrals the model already evaluates for the nitrogen +%%% profile in Canopy_Resistence_An_Evolution.m lines 55-56, with the viewing +%%% extinction Kv in place of Knit. Denominators are LAI*Fsun and LAI*Fshd +%%% exactly as defined there, so the scaling is consistent with the rest of the +%%% canopy module by construction. +%%% +%%% SCATTERING, AND WHY THE TWO-STREAM ALBEDO IS THE WRONG ONE HERE +%%% At 755 nm the leaf single-scattering albedo is high (omega_l ~ 0.87), so most +%%% interceptions scatter rather than absorb. It is tempting to reduce the +%%% extinction by sqrt(1-omega_l), which is the two-stream result (Goudriaan; +%%% Sellers 1985) -- but that applies to a DIFFUSE FLUX propagating through the +%%% medium, not to a DIRECTIONAL escape path. A scattered photon is redirected +%%% roughly isotropically, so about half of it goes back down and is lost to the +%%% sensor. Using the full albedo therefore over-credits escape. +%%% +%%% Recollision theory (Knyazikhin et al. 1998; Stenberg 2007) gives the +%%% effective albedo for escape as omega_eff = omega_l*(1 - p_recoll), where +%%% p_recoll is the probability that a scattered photon hits another leaf. For +%%% a closed canopy p_recoll is 0.5-0.7, so omega_eff is 0.26-0.44 rather than +%%% 0.87, and Kv is reduced by sqrt(1-omega_eff). +%%% +%%% This matters: with the two-stream albedo the bulk escape at LAI 3.5 comes out +%%% 0.753 and the modelled daily SIF-GPP slope is 10.27, below the observed +%%% 11.91-68.59 band from 40 FLUXNET towers (Zhang et al. 2018, Biogeosciences). +%%% With p_recoll = 0.6 the bulk escape is 0.553 and the slope is 13.59, inside +%%% the band, with r2 and RMSE essentially unchanged. Setting p_recoll = 1 +%%% recovers pure absorption and the lower bound on escape. +%%% +%%% INPUT +%%% SIF_sun, SIF_shd [W m-2 sr-1 um-1] leaf-level F755 from +%%% photosynthesis_biochemical.m +%%% LAI [-] leaf area index +%%% Kopt [-] beam extinction coefficient, from Canopy_Radiative_Transfer.m +%%% theta_v [rad] sensor zenith angle (0 = nadir; FLEX is near-nadir) +%%% omega_l [-] leaf single-scattering albedo at 755 nm, default 0.87 +%%% LADF string leaf angle distribution: 'spherical','planophile','erectophile' +%%% p_recoll [-] recollision probability, default 0.6 (closed canopy) +%%% +%%% OUTPUT +%%% fesc_sun, fesc_shd [-] escape fractions +%%% SIF_toc [W m-2 sr-1 um-1] top-of-canopy radiance +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if nargin < 5 || isempty(theta_v); theta_v = 0; end +if nargin < 6 || isempty(omega_l); omega_l = 0.87; end +if nargin < 7 || isempty(LADF); LADF = 'spherical'; end +if nargin < 8 || isempty(p_recoll); p_recoll = 0.6; end + +%%% G-function: projection of unit leaf area onto the viewing direction +switch lower(LADF) + case 'spherical'; G = 0.5; + case 'planophile'; G = cos(theta_v); + case 'erectophile'; G = 2*sin(theta_v)/pi; + otherwise; G = 0.5; +end +G = max(G,1e-3); +Kv = G./max(cos(theta_v),1e-3); +omega_eff = omega_l*(1 - min(max(p_recoll,0),1)); % recollision-corrected albedo +Kv = Kv*sqrt(max(1-omega_eff,1e-6)); % escape extinction + +if LAI <= 1e-6 + fesc_sun = 1; fesc_shd = 1; SIF_toc = 0; return +end + +Ks = max(Kopt,1e-6); +Ksv = Ks + Kv; + +Asun_esc = (1 - exp(-Ksv*LAI))/Ksv; % sunlit, escaping +Asun_tot = (1 - exp(-Ks *LAI))/Ks; % sunlit, total = LAI*Fsun +Ashd_esc = (1 - exp(-Kv *LAI))/Kv - Asun_esc; % shaded, escaping +Ashd_tot = LAI - Asun_tot; % shaded, total = LAI*Fshd + +fesc_sun = Asun_esc/max(Asun_tot,1e-9); +fesc_shd = Ashd_esc/max(Ashd_tot,1e-9); +fesc_sun = min(max(fesc_sun,0),1); +fesc_shd = min(max(fesc_shd,0),1); + +SIF_toc = SIF_sun*Asun_tot*fesc_sun + SIF_shd*Ashd_tot*fesc_shd; +end diff --git a/T&C_Code/SIF_Output.m b/T&C_Code/SIF_Output.m new file mode 100644 index 0000000..44f1c96 --- /dev/null +++ b/T&C_Code/SIF_Output.m @@ -0,0 +1,45 @@ +function [SIF_obs,mask] = SIF_Output(SIF_toc,Datam,DeltaGMT,Lon,Lat,revisit_d,overpass_h,window_h) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Subfunction SIF_Output % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Resample continuous modelled top-of-canopy SIF onto a satellite's +%%% observing schedule, so model output sits in the sensor's observation space +%%% and can be compared without further processing. +%%% +%%% Defaults are FLEX/FLORIS: 27-day repeat, ~10:00 local solar time descending +%%% node in tandem with Sentinel-3, averaged over a 1 h window. +%%% +%%% INPUT +%%% SIF_toc [W m-2 sr-1 um-1] hourly, from SIF_Escape.m +%%% Datam [Yr MO DA HR] +%%% revisit_d [d] repeat cycle, default 27 +%%% overpass_h [h] local solar time of overpass, default 10.0 +%%% window_h [h] averaging window, default 1.0 +%%% +%%% OUTPUT +%%% SIF_obs same length as SIF_toc, NaN except on overpass steps +%%% mask logical, true on overpass steps +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if nargin < 6 || isempty(revisit_d); revisit_d = 27; end +if nargin < 7 || isempty(overpass_h); overpass_h = 10.0; end +if nargin < 8 || isempty(window_h); window_h = 1.0; end + +Yr = Datam(:,1); Mo = Datam(:,2); Da = Datam(:,3); Hr = Datam(:,4); +dn = datenum(Yr,Mo,Da); +day0 = min(dn); +on_cycle = mod(dn - day0, revisit_d) == 0; + +%%% local solar time, so the comparison is at the true overpass geometry +days = [31 28 31 30 31 30 31 31 30 31 30 31]; cum = [0 cumsum(days)]; +jDay = cum(Mo)' + Da; +leap = (mod(Yr,4)==0 & (mod(Yr,100)~=0 | mod(Yr,400)==0)) & Mo>2; +jDay = jDay + double(leap); +gam = 2*pi*(jDay-1)/365; +EoT = 229.18*(0.000075 + 0.001868*cos(gam) - 0.032077*sin(gam) ... + - 0.014615*cos(2*gam) - 0.040849*sin(2*gam)); +LST = Hr + (4*(Lon - 15*DeltaGMT) + EoT)/60; + +mask = on_cycle & abs(LST - overpass_h) <= window_h/2; +SIF_obs = nan(size(SIF_toc)); +SIF_obs(mask) = SIF_toc(mask); +end diff --git a/T&C_Code/calibrate_hour_offset.m b/T&C_Code/calibrate_hour_offset.m new file mode 100644 index 0000000..780fd0b --- /dev/null +++ b/T&C_Code/calibrate_hour_offset.m @@ -0,0 +1,36 @@ +function [best, tbl] = calibrate_hour_offset(Rsw,Datam,DeltaGMT,Lon,Lat) +%%% Determine the timestamp convention from the data itself. +%%% Scans candidate offsets and returns the one minimising physically +%%% impossible hours: measurable shortwave with the sun below the horizon, or a +%%% high sun with zero shortwave. At Zurich this returns +1.0 h, reducing +%%% inconsistent hours from 6722 to 22 out of 276096. +cands = -2:0.5:2; tbl = zeros(numel(cands),3); +for k = 1:numel(cands) + C.hour_offset = cands(k); + C.f_diff=[0;0;0;0;0;0;0]; C.fvis_dir=zeros(6,1); C.fvis_dif=zeros(6,1); + C.par_dir=zeros(6,1); C.par_dif=zeros(6,1); + [~,~,~,~,~,~] = deal(0,0,0,0,0,0); + Hr = Datam(:,4) + cands(k); + cosz = local_cosz(Datam(:,1),Datam(:,2),Datam(:,3),Hr,DeltaGMT,Lon,Lat); + bad1 = sum(Rsw(:) > 20 & cosz <= 0); + bad2 = sum(cosz > 0.3 & Rsw(:) <= 0); + tbl(k,:) = [cands(k) bad1 bad2]; +end +[~,i] = min(tbl(:,2) + tbl(:,3)); +best = tbl(i,1); +end + +function cosz = local_cosz(Yr,Mo,Da,Hr,DeltaGMT,Lon,Lat) +days = [31 28 31 30 31 30 31 31 30 31 30 31]; cum = [0 cumsum(days)]; +jDay = cum(Mo)' + Da; +leap = (mod(Yr,4)==0 & (mod(Yr,100)~=0 | mod(Yr,400)==0)) & Mo>2; +jDay = jDay + double(leap); +gam = 2*pi*(jDay-1)/365; +delta = 0.006918 - 0.399912*cos(gam) + 0.070257*sin(gam) - 0.006758*cos(2*gam) ... + + 0.000907*sin(2*gam) - 0.002697*cos(3*gam) + 0.00148*sin(3*gam); +EoT = 229.18*(0.000075 + 0.001868*cos(gam) - 0.032077*sin(gam) ... + - 0.014615*cos(2*gam) - 0.040849*sin(2*gam)); +HA = ((Hr*60 + 4*(Lon - 15*DeltaGMT) + EoT)/4 - 180)*pi/180; +LatR = Lat*pi/180; +cosz = max(sin(LatR)*sin(delta) + cos(LatR)*cos(delta).*cos(HA),0); +end diff --git a/Validation/benchmark_sif_gpp.py b/Validation/benchmark_sif_gpp.py new file mode 100644 index 0000000..6fce7be --- /dev/null +++ b/Validation/benchmark_sif_gpp.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +"""benchmark_sif_gpp.py -- the chain against the published 40-tower benchmark. + +Reference: Zhang, Joiner, Alemohammad, Zhou, Gentine (2018), Biogeosciences 15, +5779-5800, "A global spatially contiguous solar-induced fluorescence (CSIF) +dataset using neural networks". Evaluating CSIF against GPP at 40 FLUXNET +tier-1 towers they report a regression slope spanning + + 11.91 to 68.59 g C m-2 day-1 per mW m-2 nm-1 sr-1 + +with per-site r2 from 0.01 to 0.93 (median 0.64) and mean RMSE 1.67 g C m-2 d-1. + +This aggregates the modelled hourly chain to those units and compares. Nothing +in the chain is fitted to the relationship. +""" +from __future__ import annotations +import argparse, os, sys +import numpy as np, scipy.io as sio +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from chain import photo_sif, escape + +LO, HI, R2_MED, RMSE_REF = 11.91, 68.59, 0.64, 1.67 + +def daily(forcing, LAI=None, Vmax=55.0, Kopt=0.5, omega=0.87, Ca=400.0, p_recoll=0.6): + d = sio.loadmat(forcing); f = lambda k: np.asarray(d[k]).ravel().astype(float) + PARB, PARD = f('PARB'), f('PARD') + Ta, ea, es, Pre = f('Ta'), f('ea'), f('esat'), f('Pre') + day = np.floor(f('D')).astype(np.int64) + PAR = PARB + PARD; Ds = np.maximum(es-ea, 0) + if LAI is None: # phenology from the site parameters + import datetime as _dt + from phenology import grass_LAI + doy = np.array([(_dt.datetime.fromordinal(int(x)-366)).timetuple().tm_yday + for x in f('D')]) + LAI, _, _, _ = grass_LAI(Ta, day, doy) + LAI = np.asarray(LAI, float)*np.ones_like(PAR) + Fsun = (1-np.exp(-Kopt*LAI))/np.maximum(Kopt*LAI, 1e-9) + PAR_sun = np.where(PAR > 0, PARB/np.maximum(Fsun,1e-6)+PARD, 0.0) + valid = (Ta > -20) & (Ta < 45); lit = (PAR > 5) & valid + Tc = np.clip(Ta, 0.1, 40) + A_s,_,F_s = photo_sif(np.where(valid, PAR_sun, 0), Ca, Tc, Ds, Pre, Vmax=Vmax) + A_h,_,F_h = photo_sif(np.where(valid, PARD, 0), Ca, Tc, Ds, Pre, Vmax=Vmax) + A_s,A_h,F_s,F_h = [np.where(lit, x, 0.0) for x in (A_s,A_h,F_s,F_h)] + om_eff = omega*(1.0 - min(max(p_recoll,0.0),1.0)) + Kv = 0.5*np.sqrt(max(1-om_eff,1e-6)); Ksv = Kopt + Kv + As = (1-np.exp(-Kopt*LAI))/Kopt + As_e = (1-np.exp(-Ksv*LAI))/Ksv + Ah = LAI - As + Ah_e = (1-np.exp(-Kv*LAI))/Kv - As_e + fs = As_e/np.maximum(As,1e-9); fh = Ah_e/np.maximum(Ah,1e-9) + _, inv = np.unique(day, return_inverse=True) + cnt = np.bincount(inv) + GPP = np.bincount(inv, weights=A_s*As + A_h*Ah)*3600*12.011e-6 # g C m-2 d-1 + SIF = np.bincount(inv, weights=F_s*As*fs + F_h*Ah*fh)/np.maximum(cnt,1) + full = cnt >= 23 + GPP, SIF = GPP[full], SIF[full] + ok = (GPP > 0.2) & (SIF > 1e-4) + sl, ic = np.linalg.lstsq(np.c_[SIF[ok], np.ones(ok.sum())], GPP[ok], rcond=None)[0] + r = np.corrcoef(SIF[ok], GPP[ok])[0,1] + rmse = float(np.sqrt(np.mean((GPP[ok]-(sl*SIF[ok]+ic))**2))) + return dict(slope=sl, r2=r**2, rmse=rmse, n=int(ok.sum()), + mGPP=GPP[ok].mean(), mSIF=SIF[ok].mean()) + +if __name__ == "__main__": + ap = argparse.ArgumentParser() + ap.add_argument("--forcing", default="Inputs/Data_Run_Zurich_Fluntern.mat") + ap.add_argument("--lai", type=float, default=None, + help="constant LAI; omit to use the site phenology") + ap.add_argument("--vmax", type=float, default=55.0) + ap.add_argument("--precoll", type=float, default=0.6) + a = ap.parse_args() + print("Benchmark: Zhang et al. 2018 BG, 40 FLUXNET tier-1 towers") + print(f" slope {LO}-{HI} g C m-2 d-1 per mW m-2 nm-1 sr-1;" + f" r2 median {R2_MED}; RMSE {RMSE_REF}\n") + R = daily(a.forcing, LAI=a.lai, Vmax=a.vmax, p_recoll=a.precoll) + print(f" modelled, n={R['n']} days, LAI={a.lai or 'site phenology'}, " + f"Vmax={a.vmax}, p_recoll={a.precoll}") + print(f" slope {R['slope']:8.2f} {'IN RANGE' if LO<=R['slope']<=HI else 'OUTSIDE'}") + print(f" r2 {R['r2']:8.4f} {'above' if R['r2']>R2_MED else 'below'} the median") + print(f" RMSE {R['rmse']:8.2f} reference {RMSE_REF}") + print(f"\n {'Vmax':>5} {'slope':>8} {'r2':>7} {'RMSE':>7} verdict") + for V in (20,30,40,55,65,80,120): + r = daily(a.forcing, LAI=a.lai, Vmax=V, p_recoll=a.precoll) + print(f" {V:5} {r['slope']:8.2f} {r['r2']:7.3f} {r['rmse']:7.2f} " + f"{'IN RANGE' if LO<=r['slope']<=HI else 'outside'}") diff --git a/Validation/chain.py b/Validation/chain.py new file mode 100644 index 0000000..1f521bd --- /dev/null +++ b/Validation/chain.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +"""chain.py -- T&C's photosynthesis + SIF chain, transcribed, driven by real forcing. + +Reproduces photosynthesis_biochemical.m (Farquhar + the Lee et al. 2015 GCB +fluorescence block) and Canopy_Resistence_An_Evolution.m line 164, then applies +SIF_Escape.m. Driven by the Zurich forcing shipped with the repo. + +The test is the emergent SIF-GPP relationship. It is not fitted anywhere in the +chain: fiF comes from the NPQ parameterisation, GPP from the Farquhar solution. +If the coupled model is right, the slope should land in the range reported from +towers and OCO-2 at 757-760 nm. +""" +from __future__ import annotations +import numpy as np + +def photo_sif(IPAR_W, Ca, Ts, Ds, Pre, Vmax=55.0, CT=3, a1=6.0, go=0.01, + rjv=1.9, Oa=210.0, Do=1000.0): + """Farquhar with the Lee et al. fluorescence block. IPAR_W in W/m2.""" + IPAR = IPAR_W*4.57 # umol photons /s/m2 + Tk = Ts + 273.15 + def arrh(c, dH): return np.exp(c - dH/(0.008314*Tk)) + Kc = arrh(38.05, 79.43); Ko = arrh(20.30, 36.38) + GAM_s = arrh(19.02, 37.83) + GAM = 0.5*np.exp(-3.3801 + 5220.0/(Tk*8.314))*Oa*Kc/Ko + kT = np.exp(26.35 - 65.33/(0.008314*Tk))/(1 + np.exp((0.71*Tk - 220.0)/(0.008314*Tk))) + Vm = Vmax*np.exp(26.35 - 65.33/(0.008314*Tk))/(1 + np.exp((0.65*Tk - 200.0)/(0.008314*Tk))) + Jmax = Vmax*rjv; Jm = Jmax*kT + FI = 0.081 # intrinsic quantum efficiency + Q = FI*IPAR # umolCO2 /s/m2 + th = 0.9 + J = (Q + Jm - np.sqrt(np.maximum((Q+Jm)**2 - 4*th*Q*Jm, 0)))/(2*th) + Rdark = 0.015*Vm + Cc = 0.7*Ca + for _ in range(30): + Cc = np.maximum(Cc, 1e-3) + JC = Vm*(Cc - GAM)/(Cc + Kc*(1 + Oa/Ko)) + JE = J*(Cc - GAM)/(4*(Cc + 2*GAM)) + A = np.minimum(JC, JE) + An = A - Rdark + gsCO2 = go + a1*An*Pre/((Cc - GAM)*(1 + Ds/Do)) + gsCO2 = np.maximum(gsCO2, go) + Cc_new = Ca - An*Pre/np.maximum(gsCO2, 1e-6) + Cc = 0.5*Cc + 0.5*np.clip(Cc_new, 1e-3, Ca) + # ---- Lee et al. 2015 fluorescence, verbatim from lines 270-291 + Jfe = A*(Cc + 2*GAM)/np.maximum(Cc - GAM, 1e-6) if CT == 3 else A + fiP0 = FI*4 + fiP = fiP0*Jfe/np.maximum(Q, 1e-9) + dls = np.clip(1 - fiP/fiP0, 0, 1) + kf = 0.05 + kd = np.maximum(0.03*Ts + 0.0773, 0.087) + kn = (6.2473*dls - 0.5944)*dls + fiF = kf/(kf + kd + kn)*(1 - fiP) + SIF = IPAR*fiF + k = 0.0375*Vmax + 8.25 + F755 = np.where(IPAR > 0, SIF/k, 0.0) + return np.maximum(A, 0.0), np.maximum(An, 0.0), np.maximum(F755, 0.0) + +def escape(LAI, Kopt=0.5, theta_v=0.0, omega=0.87, G=0.5, p_recoll=0.6): + """Escape fractions. omega_eff = omega*(1-p_recoll): a scattered photon is + redirected isotropically, so the two-stream albedo over-credits directional + escape. p_recoll = 1 recovers pure absorption.""" + omega_eff = omega*(1.0 - min(max(p_recoll, 0.0), 1.0)) + Kv = G/max(np.cos(theta_v), 1e-3)*np.sqrt(max(1-omega_eff, 1e-6)) + Ks = max(Kopt, 1e-6); Ksv = Ks + Kv + As_e = (1-np.exp(-Ksv*LAI))/Ksv; As_t = (1-np.exp(-Ks*LAI))/Ks + Ah_e = (1-np.exp(-Kv*LAI))/Kv - As_e; Ah_t = LAI - As_t + return As_e/max(As_t,1e-9), Ah_e/max(Ah_t,1e-9), As_t, Ah_t diff --git a/Validation/phenology.py b/Validation/phenology.py new file mode 100644 index 0000000..bc84dff --- /dev/null +++ b/Validation/phenology.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python3 +"""phenology.py -- prescribed grass LAI from the site's own T&C parameters. + +MOD_PARAM_ZURICH_SMA.m gives, for the low-vegetation layer: + aSE_L = 2 grass species + Tlo_L = 0.0 mean temperature for leaf onset [C] + Tls_L = NaN no temperature-driven leaf shed + dmg_L = 20 days of maximum growth + LAI_min_L = 0.1 minimum LAI + Sl_L = 0.035 specific leaf area [m2/gC] + +Holding LAI at 4 year-round contradicts LAI_min_L = 0.1 by a factor of forty in +dormancy, and assigns a full canopy to grass that is not there. In the SIF-GPP +regression that inflates GPP at the low-SIF end and flattens the slope, which is +the wrong direction for comparison against tower data. + +This drives LAI from a running-mean temperature threshold at Tlo, with a +logistic build over dmg days and autumn senescence, bounded by LAI_min and a +prescribed peak. It is a prescribed phenology, not T&C's prognostic +VEGETATION_DYNAMIC, and it is used here only so the benchmark is run on a +canopy consistent with the site's own parameters. +""" +from __future__ import annotations +import numpy as np + +def running_mean_daily(Ta_hourly, day_idx, window=7): + u, inv = np.unique(day_idx, return_inverse=True) + cnt = np.bincount(inv) + Td = np.bincount(inv, weights=Ta_hourly)/np.maximum(cnt, 1) + k = np.ones(window)/window + Tsm = np.convolve(np.r_[np.repeat(Td[0], window), Td], k, mode="same")[window:] + return u, Td, Tsm[:len(Td)] + +def grass_LAI(Ta_hourly, day_idx, doy_hourly, Tlo=0.0, dmg=20, LAI_min=0.1, + LAI_max=3.5, sen_doy=270, sen_len=45): + """Daily LAI, mapped back onto the hourly index.""" + u, Td, Tsm = running_mean_daily(Ta_hourly, day_idx) + _, inv = np.unique(day_idx, return_inverse=True) + doy_d = np.bincount(inv, weights=doy_hourly)/np.maximum(np.bincount(inv), 1) + LAI_d = np.full(len(u), LAI_min) + grow = 0.0 + for i in range(len(u)): + active = (Tsm[i] > Tlo) and (doy_d[i] < sen_doy) + if active: + grow = min(grow + 1.0/dmg, 1.0) + else: + grow = max(grow - 1.0/sen_len, 0.0) + LAI_d[i] = LAI_min + (LAI_max - LAI_min)*grow + return LAI_d[inv], LAI_d, u, doy_d diff --git a/Validation/run_dehai.py b/Validation/run_dehai.py new file mode 100644 index 0000000..3d335ed --- /dev/null +++ b/Validation/run_dehai.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python3 +"""run_dehai.py -- validate the chain at DE-Hai (Hainich) against measurements. + +WHY THIS SITE. The Zurich benchmark compares against MODELLED bands: the four +shortwave components shipped in the repository are 98.7% determined by (cloud +fraction, cos z), so they are a partitioning model's output, not instrument +readings. DE-Hai measures diffuse shortwave (SW_DIF) directly and reports +eddy-covariance GPP, so two links in the chain become testable against +observation rather than against another model: + + 1 the direct/diffuse partition, against measured SW_DIF + 2 simulated GPP, against GPP_NT_VUT_REF and GPP_DT_VUT_REF + +SIF itself remains unvalidated here -- DE-Hai carries no fluorescence +spectrometer -- but it is then driven by a radiation partition and a +photosynthesis rate that have each been checked against measurement, which is a +different claim from the aggregate consistency the 40-tower benchmark gives. + +DATA. ICOS Carbon Portal, ICOSETC_DE-Hai_FLUXNET_FLUXMET_HH_*.csv, or the +FLUXNET2015 FULLSET product. Half-hourly; this aggregates to hourly. + + python run_dehai.py --csv ICOSETC_DE-Hai_FLUXNET_FLUXMET_HH_2000-2025_v1.3.csv + +SITE. Hainich, Germany. Deciduous broadleaf (Fagus sylvatica), 51.0792 N, +10.4522 E, 430 m, UTC+1. Peak LAI 5-6, Vcmax25 for European beech 45-70. +""" +from __future__ import annotations +import argparse, os, sys, warnings +warnings.filterwarnings("ignore", category=RuntimeWarning) +import numpy as np +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), + "..", "Forcing_Prep")) +from chain import photo_sif +from partition import solar_geometry, erbs_diffuse_fraction, SOLAR_CONST + +LAT, LON, ELEV, GMT = 51.0792, 10.4522, 430.0, 1.0 +MISSING = -9999.0 + +def read_icos(path): + raw = np.genfromtxt(path, delimiter=",", names=True, dtype=None, encoding="utf-8") + c = raw.dtype.names + def col(*names): + for n in names: + if n in c: + v = np.asarray(raw[n], float) + return np.where(v <= MISSING + 1, np.nan, v) + return None + ts = np.asarray(raw["TIMESTAMP_START"], str) + return dict( + yr=np.array([int(s[0:4]) for s in ts]), mo=np.array([int(s[4:6]) for s in ts]), + dy=np.array([int(s[6:8]) for s in ts]), + hr=np.array([int(s[8:10]) + int(s[10:12])/60 for s in ts], float), + SW_IN=col("SW_IN_F","SW_IN"), SW_DIF=col("SW_DIF"), + PPFD_IN=col("PPFD_IN"), PPFD_DIF=col("PPFD_DIF"), + TA=col("TA_F","TA"), VPD=col("VPD_F","VPD"), PA=col("PA_F","PA"), + GPP_NT=col("GPP_NT_VUT_REF"), GPP_DT=col("GPP_DT_VUT_REF"), + LAI=col("LAI")) + +def calibrate_offset(SW, yr, mo, dy, hr): + best = None + for off in np.arange(-2.0, 2.01, 0.5): + h, _ = solar_geometry(yr, mo, dy, hr+off, LAT, LON, GMT) + cz = np.maximum(np.sin(h), 0.0) + bad = int((np.nan_to_num(SW) > 20).__and__(cz <= 0).sum() + + ((cz > 0.3) & (np.nan_to_num(SW) <= 0)).sum()) + if best is None or bad < best[1]: best = (off, bad) + return best + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--csv", required=True) + ap.add_argument("--vmax", type=float, default=55.0) + ap.add_argument("--laimax", type=float, default=5.5) + ap.add_argument("--precoll", type=float, default=0.6) + ap.add_argument("--kopt", type=float, default=0.5) + a = ap.parse_args() + D = read_icos(a.csv) + n = len(D["yr"]); print(f"{n} records, {D['yr'].min()}-{D['yr'].max()}") + for k in ("SW_DIF","GPP_NT","PPFD_IN"): + print(f" {k:8} {'present' if D[k] is not None else 'ABSENT'}" + + (f", {np.isfinite(D[k]).mean()*100:.1f}% finite" if D[k] is not None else "")) + + off, bad = calibrate_offset(D["SW_IN"], D["yr"], D["mo"], D["dy"], D["hr"]) + print(f"\n timestamp offset {off:+.1f} h ({bad} physically impossible hours)") + h, r = solar_geometry(D["yr"], D["mo"], D["dy"], D["hr"]+off, LAT, LON, GMT) + cz = np.maximum(np.sin(h), 0.0); I0 = SOLAR_CONST*r*cz + SW = D["SW_IN"]; kt = np.where(I0 > 1, SW/np.maximum(I0, 1e-9), 0.0) + + # ---- TEST 1: direct/diffuse partition against MEASURED SW_DIF + if D["SW_DIF"] is not None: + obs_fd = D["SW_DIF"]/np.maximum(SW, 1e-9) + m = np.isfinite(obs_fd) & (SW > 20) & (obs_fd >= 0) & (obs_fd <= 1) + fd_erbs = erbs_diffuse_fraction(kt) + rmse = lambda p: float(np.sqrt(np.nanmean((p[m]-obs_fd[m])**2))) + print(f"\n TEST 1 diffuse fraction vs MEASURED SW_DIF, n={m.sum()}") + print(f" Erbs (1982) RMSE {rmse(fd_erbs):.4f} bias " + f"{np.nanmean((fd_erbs-obs_fd)[m]):+.4f}") + from tcpart import TCPartition, predictors, _design_fd, fit_logit + ktp, mm, lm, czp, pers = predictors(SW, h, r, np.exp(-ELEV/8434.0)) + idx = np.arange(n) # split by position, not by year: + cut = idx[m][int(0.6*m.sum())] if m.sum() > 10 else n # a year-based split + tr = m & (idx <= cut) # is empty for single-year files + b = fit_logit(_design_fd(ktp[tr], lm[tr], czp[tr], pers[tr]), obs_fd[tr], SW[tr]) + fd_fit = 1/(1+np.exp(-np.clip(_design_fd(ktp,lm,czp,pers)@b,-40,40))) + te = m & ~tr + print(f" calibrated, held out RMSE " + f"{float(np.sqrt(np.nanmean((fd_fit[te]-obs_fd[te])**2))):.4f} " + f"bias {np.nanmean((fd_fit-obs_fd)[te]):+.4f} n={te.sum()}") + fd_use = fd_fit + else: + print("\n TEST 1 skipped: no SW_DIF column") + fd_use = erbs_diffuse_fraction(kt) + + # ---- LAI: measured if present, else beech phenology + if D["LAI"] is not None and np.isfinite(D["LAI"]).mean() > 0.3: + LAI = np.where(np.isfinite(D["LAI"]), D["LAI"], 0.5); src = "measured" + else: + doy = np.array([sum([31,28,31,30,31,30,31,31,30,31,30,31][:mo-1])+dd + for mo, dd in zip(D["mo"], D["dy"])], float) + g = np.clip((doy-110)/30, 0, 1)*np.clip((300-doy)/30, 0, 1) + LAI = 0.3 + (a.laimax-0.3)*g; src = "beech phenology, DOY 110-300" + print(f"\n LAI: {src}, range {np.nanmin(LAI):.2f}-{np.nanmax(LAI):.2f}") + + # ---- TEST 2: simulated GPP against tower GPP + PAR = 0.46*SW + PARD = PAR*fd_use; PARB = PAR - PARD + Fsun = (1-np.exp(-a.kopt*LAI))/np.maximum(a.kopt*LAI, 1e-9) + Ta = D["TA"]; Ds = D["VPD"]*100.0 if D["VPD"] is not None else 1000.0 + Pre = D["PA"]*10.0 if D["PA"] is not None else 1013.0 + ok = np.isfinite(Ta) & np.isfinite(SW) + Tc = np.clip(np.nan_to_num(Ta, nan=10.0), 0.1, 40) + Dsc = np.clip(np.nan_to_num(Ds, nan=1000.0), 10, 6000) + Prc = np.nan_to_num(Pre, nan=1013.0) + lit = ok & (PAR > 5) + A_s,_,F_s = photo_sif(np.where(lit, PARB/np.maximum(Fsun,1e-6)+PARD, 0), 400., Tc, Dsc, Prc, Vmax=a.vmax) + A_h,_,F_h = photo_sif(np.where(lit, PARD, 0), 400., Tc, Dsc, Prc, Vmax=a.vmax) + A_s,A_h,F_s,F_h = [np.where(lit, x, 0.0) for x in (A_s,A_h,F_s,F_h)] + om = 0.87*(1.0 - min(max(a.precoll,0),1)); Kv = 0.5*np.sqrt(max(1-om,1e-6)); Ksv = a.kopt+Kv + As = (1-np.exp(-a.kopt*LAI))/a.kopt; As_e = (1-np.exp(-Ksv*LAI))/Ksv + Ah = LAI-As; Ah_e = (1-np.exp(-Kv*LAI))/Kv - As_e + fs = As_e/np.maximum(As,1e-9); fh = Ah_e/np.maximum(Ah,1e-9) + GPP_mod = A_s*As + A_h*Ah + SIF_toc = F_s*As*fs + F_h*Ah*fh + synth = os.path.basename(a.csv).lower().find("synth") >= 0 + if synth: + print("\n NOTE: filename contains 'synth'. TEST 2 compares against generated") + print(" GPP and is a harness check only; the numbers carry no meaning.") + for nm in ("GPP_NT","GPP_DT"): + if D[nm] is None: continue + g = D[nm]; mm2 = lit & np.isfinite(g) & (g > -5) + if mm2.sum() < 100: continue + bias = float(np.nanmean((GPP_mod-g)[mm2])) + rm = float(np.sqrt(np.nanmean((GPP_mod-g)[mm2]**2))) + rr = float(np.corrcoef(GPP_mod[mm2], g[mm2])[0,1]) + sl = np.linalg.lstsq(np.c_[g[mm2], np.ones(mm2.sum())], GPP_mod[mm2], rcond=None)[0][0] + print(f"\n TEST 2 simulated GPP vs tower {nm}, n={mm2.sum()}") + print(f" r2 {rr**2:.4f} slope {sl:.3f} bias {bias:+.3f} RMSE {rm:.3f} umol m-2 s-1") + print(f" mean modelled {GPP_mod[mm2].mean():.2f} mean tower {g[mm2].mean():.2f}") + + # ---- daily SIF-GPP, in the 40-tower benchmark units + day = D["yr"]*10000 + D["mo"]*100 + D["dy"] + _, inv = np.unique(day, return_inverse=True); cnt = np.bincount(inv) + dt_h = 0.5 if np.median(np.diff(D["hr"][:48])) < 0.75 else 1.0 + G = np.bincount(inv, weights=np.nan_to_num(GPP_mod))*3600*dt_h*12.011e-6 + S = np.bincount(inv, weights=np.nan_to_num(SIF_toc))/np.maximum(cnt,1) + full = cnt >= (23/dt_h); G, S = G[full], S[full] + o = (G > 0.2) & (S > 1e-4) + if o.sum() > 100: + sl, ic = np.linalg.lstsq(np.c_[S[o], np.ones(o.sum())], G[o], rcond=None)[0] + rr = np.corrcoef(S[o], G[o])[0,1] + print(f"\n daily SIF-GPP, n={o.sum()} days") + print(f" slope {sl:.2f} r2 {rr**2:.4f} " + f"{'IN RANGE' if 11.91<=sl<=68.59 else 'outside'} (11.91-68.59)") + +if __name__ == "__main__": + main() diff --git a/Validation/sif_gpp.py b/Validation/sif_gpp.py new file mode 100644 index 0000000..be4b454 --- /dev/null +++ b/Validation/sif_gpp.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +"""sif_gpp.py -- the emergent SIF-GPP relationship of the T&C chain. + +Runs photosynthesis_biochemical.m's Farquhar + Lee et al. 2015 fluorescence +block, the sunlit/shaded scaling of Canopy_Resistence_An_Evolution.m line 164, +and SIF_Escape.m, on real forcing. Nothing in the chain is fitted to the +SIF-GPP relationship, so the relationship is a prediction. + + python sif_gpp.py [--forcing PATH] [--lai 4] [--vmax 55] +""" +from __future__ import annotations +import argparse, sys, os +import numpy as np, scipy.io as sio, datetime as dt +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from chain import photo_sif, escape + +def load(path): + d = sio.loadmat(path); f = lambda k: np.asarray(d[k]).ravel().astype(float) + D = f('D'); n = len(D); mo = np.empty(n, int) + for i, dv in enumerate(D): + mo[i] = (dt.datetime.fromordinal(int(dv)-366)+dt.timedelta(days=float(dv)%1)).month + return dict(PARB=f('PARB'), PARD=f('PARD'), Ta=f('Ta'), ea=f('ea'), + esat=f('esat'), Pre=f('Pre'), mo=mo) + +def run(F, LAI=4.0, Vmax=55.0, Kopt=0.5, omega=0.87, Ca=400.0, theta_v=0.0): + PAR = F['PARB']+F['PARD']; Ds = np.maximum(F['esat']-F['ea'], 0) + m = (PAR > 5) & (F['Ta'] > 0) & (F['Ta'] < 40) + Fsun = (1-np.exp(-Kopt*LAI))/(Kopt*LAI) + PAR_sun = F['PARB'][m]/max(Fsun,1e-6) + F['PARD'][m] + PAR_shd = F['PARD'][m] + A_s,_,Fl_s = photo_sif(PAR_sun, Ca, F['Ta'][m], Ds[m], F['Pre'][m], Vmax=Vmax) + A_h,_,Fl_h = photo_sif(PAR_shd, Ca, F['Ta'][m], Ds[m], F['Pre'][m], Vmax=Vmax) + fs, fh, As, Ah = escape(LAI, Kopt, theta_v, omega) + GPP = A_s*As + A_h*Ah + SIF = Fl_s*As*fs + Fl_h*Ah*fh + ok = (GPP > 0.5) & (SIF > 1e-5) + sl, ic = np.linalg.lstsq(np.c_[SIF[ok], np.ones(ok.sum())], GPP[ok], rcond=None)[0] + r = np.corrcoef(SIF[ok], GPP[ok])[0,1] + return dict(slope=sl, r2=r**2, n=int(ok.sum()), fesc_sun=fs, fesc_shd=fh, + fesc_bulk=(As*fs+Ah*fh)/(As+Ah), SIF=SIF, GPP=GPP, ok=ok, + mo=F['mo'][m], mean_SIF=SIF[ok].mean(), mean_GPP=GPP[ok].mean()) + +if __name__ == "__main__": + ap = argparse.ArgumentParser() + ap.add_argument("--forcing", default="/home/claude/TeC/Inputs/Data_Run_Zurich_Fluntern.mat") + ap.add_argument("--lai", type=float, default=4.0) + ap.add_argument("--vmax", type=float, default=55.0) + a = ap.parse_args() + F = load(a.forcing) + R = run(F, LAI=a.lai, Vmax=a.vmax) + print(f"n = {R['n']} LAI = {a.lai} Vmax = {a.vmax}") + print(f" escape: sunlit {R['fesc_sun']:.4f} shaded {R['fesc_shd']:.4f} bulk {R['fesc_bulk']:.4f}") + print(f" GPP = {R['slope']:.2f} x SIF r2 = {R['r2']:.4f}") + print(f" mean SIF {R['mean_SIF']:.3f} W m-2 sr-1 um-1 mean GPP {R['mean_GPP']:.2f} umol m-2 s-1") + print("\n SENSITIVITY") + print(f" {'LAI':>5} {'slope':>7} {'r2':>7} {'bulk fesc':>10}") + for L in (1,2,3,4,6,8): + r = run(F, LAI=L, Vmax=a.vmax) + print(f" {L:5} {r['slope']:7.2f} {r['r2']:7.3f} {r['fesc_bulk']:10.3f}") + print(f" {'Vmax':>5} {'slope':>7} {'r2':>7}") + for V in (30,40,55,80,120): + r = run(F, LAI=a.lai, Vmax=V) + print(f" {V:5} {r['slope']:7.2f} {r['r2']:7.3f}")