diff --git a/README.md b/README.md index 3dae7b4c..f4fcdc66 100644 --- 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 index 2168b929..e656b2f4 100755 --- 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!")