Merge pull request #95 from angehung5/main

Add pavd and daily LAI
noaa-oar-arl · Oct 19, 2023 · 1caad6f · 1caad6f
2 parents 1146880 + f7856c8
commit 1caad6f
Show file tree

Hide file tree

Showing 2 changed files with 88 additions and 23 deletions.
diff --git a/README.md b/README.md
@@ -150,10 +150,11 @@ The Canopy-App input data in [Table 2](#table-2-canopy-app-required-input-variab
 | `hpbl`                           | Height of the planetary boundary layer (m)  | UFS NOAA/GFSv16                                    |
 | `prate_ave`                      | Average mass precipitation rate (kg m-2 s-1) | UFS NOAA/GFSv16                                   |
 | **External Canopy Variables**    | **Variable Description and Units**          | **Data Source/Reference (if necessary)**           |
-| `fh`                             | Forest canopy height (m)                    | Fused GEDI/Landsat data. Data Period=2020. ([Potapov et al., 2020](https://doi.org/10.1016/j.rse.2020.112165)) |
-| `clu`                            | Forest clumping index (dimensionless)       | GriddingMachine/MODIS. Data Period=2001-2017 Climatology. ([Wei et al., 2019](https://doi.org/10.1016/j.rse.2019.111296)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
-| `lai`                            | Leaf area index (m2/m2)                     | VIIRS-NPP. Data Period=2018-2020 Climatology. ([Myneni 2018](https://doi.org/10.5067/VIIRS/VNP15A2H.001)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
-| `ffrac`                          | Forest fraction (dimensionless)             | Based on [VIIRS GVF-NPP](https://www.star.nesdis.noaa.gov/jpss/gvf.php) and GriddingMachine/MODIS FFRAC/FVC from Terra.  Data Period=2020. ([DiMiceli et al., 2022](https://doi.org/10.5067/MODIS/MOD44B.061)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
+| `fh`                             | Forest canopy height (m)                    | Fused GEDI/Landsat data. Data Period=2020. Data frequency=Annual. ([Potapov et al., 2020](https://doi.org/10.1016/j.rse.2020.112165)) |
+| `clu`                            | Forest clumping index (dimensionless)       | GriddingMachine/MODIS. Data Period=2001-2017 Climatology. Data frequency=Monthly. ([Wei et al., 2019](https://doi.org/10.1016/j.rse.2019.111296)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
+| `lai`                            | Leaf area index (m2/m2)                     | VIIRS-NPP. Data Period=2020. Data frequency=Daily, interpolated from original 8-day product. ([Myneni 2018](https://doi.org/10.5067/VIIRS/VNP15A2H.001)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
+| `ffrac`                          | Forest fraction (dimensionless)             | Based on [VIIRS GVF-NPP](https://www.star.nesdis.noaa.gov/jpss/gvf.php) and GriddingMachine/MODIS FFRAC/FVC from Terra.  Data Period=2020. Data frequency=Annual. ([DiMiceli et al., 2022](https://doi.org/10.5067/MODIS/MOD44B.061)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
+| `pavd`                           | Plant area volume density (m2/m3)           | [GEDI product from North Arizona University](https://goetzlab.rc.nau.edu/index.php/gedi/). Data Period=201904-202212 Climatology. Data frequency=Annual. Three dimensional structure of plant area volume density with 14 vertical layers from the surface (0 m) to 70 m above ground level. Data at each layer represents the average pavd within certain height range (e.g. 0 - 5 m for first layer). |
 | **Other External Variables**     | **Variable Description and Units**          | **Data Source/Reference (if necessary)**           |
 | `frp`                            | Total Fire Radiative Power (MW/grid cell area) | [NOAA/NESDIS GBBEPx](https://www.ospo.noaa.gov/Products/land/gbbepx/) |
 | `csz`                            | Cosine of the solar zenith angle (dimensionless) | [Based on Python Pysolar](https://pysolar.readthedocs.io/en/latest/) |

diff --git a/python/global_data_process.py b/python/global_data_process.py
@@ -1,5 +1,6 @@
 """
-Created on Sun Jun 4 15:45:29 2023
+Created on Sun Jun 4 2023
+Updated on Mon Oct 16 2023: Use daily gfs.canopy files
 
 Author: Wei-Ting Hung
 """
@@ -65,7 +66,7 @@
 f_met = (
     path + "/gfs.t" + HHI + "z." + YY + MM + DD + ".sfcf0" + HH + ".nc"
 )  # gfs met file
-f_can = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc"  # canopy file
+f_can = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".sfcf000.nc"  # canopy file
 f_output = (
     path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".global.f0" + HH + ".nc"
 )  # output file
@@ -88,7 +89,7 @@
             + ".nc "
         )
 elif frp_src == 2:  # climatological frp
-    f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc"
+    f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".sfcf000.nc"
 
 
 # required variables
@@ -109,7 +110,7 @@
     "hpbl",
     "prate_ave",
 ]
-canlist = ["lai", "clu", "ffrac", "fh", "mol", "csz", "frp", "href"]
+canlist = ["lai", "clu", "ffrac", "fh", "pavd", "mol", "csz", "frp", "href"]
 
 
 # constants
@@ -149,12 +150,35 @@ def mapping(xgrid, ygrid, data, xdata, ydata, map_method, fvalue):
 def read_gfs_climatology(filename, lat, lon, varname):
     readin = Dataset(filename)
 
-    # map to met grids
-    yt = readin["lat"][:]
-    xt = readin["lon"][:]
-    data = np.squeeze(readin[varname][0, :, :])
+    if varname == "pavd":
+        # map to met grids
+        yt = readin["lat"][:]
+        xt = readin["lon"][:]
+        data = np.squeeze(readin[varname][:])
+
+        DATA = np.empty([data.shape[0], lat.shape[0], lat.shape[1]])
+
+        for ll in np.arange(data.shape[0]):
+            DATA[ll, :, :] = mapping(
+                lat,
+                lon,
+                data[ll, :, :].flatten(),
+                yt.flatten(),
+                xt.flatten(),
+                "linear",
+                np.nan,
+            )
+
+    else:
+        # map to met grids
+        yt = readin["lat"][:]
+        xt = readin["lon"][:]
+        data = np.squeeze(readin[varname][0, :, :])
+
+        DATA = mapping(
+            lat, lon, data.flatten(), yt.flatten(), xt.flatten(), "linear", np.nan
+        )
 
-    DATA = mapping(lat, lon, data.flatten(), yt.flatten(), xt.flatten(), "linear", np.nan)
     DATA[np.isnan(DATA)] = 0
     DATA[DATA < 0] = 0
     return DATA
@@ -226,7 +250,8 @@ def read_frp_local(filename, lat, lon, fill_value):
             + "z."
             + YY
             + MM
-            + "01.sfcf000.nc",
+            + DD
+            + ".sfcf000.nc",
         ]
     )
     if os.path.isfile(f_can) is True:
@@ -262,7 +287,7 @@ def read_frp_local(filename, lat, lon, fill_value):
     else:
         print("---- No available FRP file. Switch to Climatology FRP...")
         frp_src = 2
-        f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc"
+        f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".sfcf000.nc"
 
 if frp_src == 2:  # 12 month climatology frp
     if os.path.isfile(f_frp) is True:
@@ -281,7 +306,8 @@ def read_frp_local(filename, lat, lon, fill_value):
                 + "z."
                 + YY
                 + MM
-                + "01.sfcf000.nc",
+                + DD
+                + ".sfcf000.nc",
             ]
         )
         if os.path.isfile(f_met) is True:
@@ -354,6 +380,11 @@ def read_frp_local(filename, lat, lon, fill_value):
         ATT = ["Canopy height above the surface", "m", fill_value]
         DATA = read_gfs_climatology(f_can, lat, lon, "fh")
 
+    elif varname == "pavd":
+        ATTNAME = ["long_name", "units", "missing_value"]
+        ATT = ["Plant area volume density profile", "m2/m3", fill_value]
+        DATA = read_gfs_climatology(f_can, lat, lon, "pavd")
+
     elif varname == "mol":
         # Reference:
         # Essa 1999, ESTIMATION OF MONIN-OBUKHOV LENGTH USING RICHARDSON AND BULK RICHARDSON NUMBER
@@ -409,15 +440,48 @@ def read_frp_local(filename, lat, lon, fill_value):
     print(ATT)
 
     # adding to output file
-    output = Dataset(f_output, "a")
+    if varname == "pavd":
+        output = Dataset(f_output, "a")
+        output.createDimension("level", DATA.shape[0])
+
+        var_bot = output.createVariable("layer_bottom", "i4", ("level",))
+        var_top = output.createVariable("layer_top", "i4", ("level",))
+        var = output.createVariable(
+            varname,
+            "float",
+            ("time", "level", "grid_yt", "grid_xt"),
+            fill_value=fill_value,
+        )
 
-    var = output.createVariable(
-        varname, "float", ("time", "grid_yt", "grid_xt"), fill_value=fill_value
-    )
-    write_varatt(var, ATTNAME, ATT)
-    var[:] = DATA
+        write_varatt(var, ATTNAME, ATT)
+        write_varatt(
+            var_bot,
+            ["long_name", "units"],
+            ["height of the layer bottom above the ground", "m"],
+        )
+        write_varatt(
+            var_top,
+            ["long_name", "units"],
+            ["height of the layer top above the ground", "m"],
+        )
+
+        var_bot[:] = np.arange(0, 65 + 1, 5)
+        var_top[:] = np.arange(5, 70 + 1, 5)
+        var[:] = DATA
+
+        output.close()
+        del [var_bot, var_top]
+
+    else:
+        output = Dataset(f_output, "a")
+
+        var = output.createVariable(
+            varname, "float", ("time", "grid_yt", "grid_xt"), fill_value=fill_value
+        )
+        write_varatt(var, ATTNAME, ATT)
+        var[:] = DATA
 
-    output.close()
+        output.close()
 
     print("---- " + varname + " complete!")