From 21b2f583b2080aaadd10f7afe0c5f93e6116d886 Mon Sep 17 00:00:00 2001
From: "Woude, Auke van der" <auke.vanderwoude@wur.nl>
Date: Wed, 17 Aug 2022 10:23:16 +0200
Subject: [PATCH] Write TER and GPP, as well as NEE

---
 da/preprocessing/upscale_bio.py | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/da/preprocessing/upscale_bio.py b/da/preprocessing/upscale_bio.py
index 3316f1c8..7a3d7190 100644
--- a/da/preprocessing/upscale_bio.py
+++ b/da/preprocessing/upscale_bio.py
@@ -33,11 +33,15 @@ class SiBUpscaler(Regional):
         self.sibfile = self.sibdir + sibfile
         self.sibfile_avg = self.sibfile.replace('lu.qp3', 'g.qp2')
         logging.debug('Reading SiB data from {}'.format(self.sibfile))
-        nee = self.load_biosphere_emissions() 
-        outfile = self.write(nee)
+        datadict = {}
+        for v in ('nee', 'resp_tot', 'assim'):
+            logging.info(f'Getting high-res {v}')
+            data = self.load_biosphere_emissions(v)
+            datadict[v] = data
+        outfile = self.write(datadict)
         return outfile, self.name
 
-    def load_biosphere_emissions(self):
+    def load_biosphere_emissions(self, fluxname):
         """Function that reads in the biosphere emissions and clips them to the right format"""
     
         # Read in the high-res pft data
@@ -49,7 +53,7 @@ class SiBUpscaler(Regional):
     
         # Load the raw SiBdata
         with nc.Dataset(self.sibfile) as ds:
-            nee = ds['nee']
+            nee = ds[fluxname]
             self.sibunits = nee.units
             sibtimes = nc.num2date(ds['time'][:], ds['time'].units) - dtm.timedelta(hours=0.5)
             startindex = np.argmin(abs(sibtimes - self.times[0]))
@@ -85,7 +89,11 @@ class SiBUpscaler(Regional):
         # Also load the pft-averaged NEE, for filling:
         logging.debug('Loading pft-averaged NEE')
         with nc.Dataset(self.sibfile_avg) as ds:
-            nee_averaged = ds['nee']
+            if fluxname == 'resp_tot':
+                fluxname = 'reco'
+            elif fluxname == 'assim':
+                fluxname = 'gpp'
+            nee_averaged = ds[fluxname]
             nee_averaged = nee_averaged[startindex:stopindex, :].astype(np.float32)
             nee_averaged = np.flip(st.convert_2d(nee_averaged, latsib, lonsib), axis=-2)
             nee_averaged = nee_averaged[:,
-- 
GitLab