change MODIS TVF to VIIRS GVF

README.md

@@ -162,7 +162,7 @@ The Canopy-App input data in [Table 2](#table-2-canopy-app-required-input-variab
 | `ch`                             | Canopy height (m)                    | Globally extended GEDI data. Data Period=2020. Data frequency=Annual. ([Lang et al., 2023](https://doi.org/10.1038/s41559-023-02206-6)) |
 | `clu`                            | Canopy 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). |
-| `canfrac`                          | Canopy fraction (dimensionless)             | Based on [MODIS VCF](https://doi.org/10.5067/MODIS/MOD44B.061). Data Period=2020. Data frequency=Annual. Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
+| `canfrac`                          | Canopy green vegetation fraction (dimensionless)             | Based on [VIIRS GVF](https://www.star.nesdis.noaa.gov/jpss/gvf.php). Data Period=2020. Data frequency=Monthly. 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). |
 | `lev`                            | Height AGL for levels associated with optional pavd (or other canopy profile) inputs (m)                                  | Same as for GEDI PAVD (or other canopy profile inputs) above                                    |
 | **Other External Variables**     | **Variable Description and Units**          | **Data Source/Reference (if necessary)**           |
@@ -188,12 +188,14 @@ Hourly gridded GFSv16 data is available from March 23, 2021 - Current Day and is
 /groups/ESS/whung/canopy_wind/gfsv16_test_data/test_2022
 ```
 
+
 **For NOAA Hera users, daily global canopy files for 2022 at 12 UTC are available at**
 
 ```
 /scratch1/RDARCH/rda-arl-gpu/Wei-ting.Hung/Global_canopy/canopy_app_2022
 ```
 
+
 **For NOAA HPSS users (e.g., Hera or WCOSS2), hourly operational GFSv16 meteorology files are archived at**
 
 ```

input/point_file_20220630.sfcf023.txt

@@ -1,2 +1,2 @@
 lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave,soilw1,soilw2,soilw3,soilw4,wilt
-34.03,272.11,20.8692,-0.7485,0.3286,0.5214,3.7059,4,0.8815,0.1107,0.0416,0.8260,12.2366,0.0000,10.00,3,99731.6094,10.1607,-9.8159,292.9766,293.2202,0.0146,14.0288,0.0001,0.2591,0.2053,0.2134,0.2132,0.0470
+34.03,272.11,20.8692,-0.7485,0.3286,0.5214,3.7059,4,0.2966,0.1107,0.0416,0.8260,12.2366,0.0000,10.00,3,99731.6094,10.1607,-9.8159,292.9766,293.2202,0.0146,14.0288,0.0001,0.2591,0.2053,0.2134,0.2132,0.0470

input/point_file_20220701.sfcf000.txt

@@ -1,2 +1,2 @@
 lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave,soilw1,soilw2,soilw3,soilw4,wilt
-34.03,272.11,20.8692,0.1580,0.8375,0.5105,3.6504,4,0.8815,0.1548,0.2357,0.8260,-181.1201,0.0000,10.00,3,99809.6719,122.7623,1.8252,294.9718,294.9180,0.0157,85.5025,0.0000,0.2487,0.2256,0.2283,0.2307,0.0470
+34.03,272.11,20.8692,0.1580,0.8375,0.5105,3.6504,4,0.2966,0.1548,0.2357,0.8260,-181.1201,0.0000,10.00,3,99809.6719,122.7623,1.8252,294.9718,294.9180,0.0157,85.5025,0.0000,0.2487,0.2256,0.2283,0.2307,0.0470

input/point_file_20220701.sfcf001.txt

@@ -1,2 +1,2 @@
 lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave,soilw1,soilw2,soilw3,soilw4,wilt
-34.03,272.11,20.8692,-0.0656,0.8462,0.5105,3.6504,4,0.8815,0.1866,0.4319,0.8260,-16.3281,0.0000,10.00,3,99854.1172,336.8927,35.6389,296.9243,296.5100,0.0162,234.9751,0.0000,0.2477,0.2259,0.2283,0.2307,0.0470
+34.03,272.11,20.8692,-0.0656,0.8462,0.5105,3.6504,4,0.2966,0.1866,0.4319,0.8260,-16.3281,0.0000,10.00,3,99854.1172,336.8927,35.6389,296.9243,296.5100,0.0162,234.9751,0.0000,0.2477,0.2259,0.2283,0.2307,0.0470

python/environment.yml

@@ -22,4 +22,4 @@ dependencies:
   # Extra
   - ipython
   - pysolar
-  - scipy
+  - monet

python/global_data_process.py

@@ -3,6 +3,7 @@
 Updated on Mon Oct 16 2023: Use daily gfs.canopy files
 Updated on Tue Nov 7  2023: Enable multiple times as user argument
 Updated on Tue Apr 2  2024: Remove wget functions, all data must be from local files
+Updated on Fri May 31 2024: Replace scipy griddata with monet (pyresample)
 
 Author: Wei-Ting Hung
 """
@@ -11,10 +12,11 @@
 import sys
 from datetime import datetime, timedelta, timezone
 
+import monet
 import numpy as np
+import xarray as xr
 from netCDF4 import Dataset
 from pysolar.solar import get_altitude
-from scipy.interpolate import griddata
 
 """User Arguments"""
 # Time: yyyymmddhhfff
@@ -103,65 +105,44 @@ def write_varatt(var, attname, att):
             var.setncattr(attname[X], att[X])
 
 
-def mapping(xgrid, ygrid, data, xdata, ydata, map_method, fvalue):
-    output = griddata(
-        (xdata, ydata), data, (xgrid, ygrid), method=map_method, fill_value=fvalue
+def read_gfs_climatology(filename, basefile, varname):
+    readin = xr.open_dataset(filename)
+    readin = readin.set_coords(["lat", "lon"]).rename(
+        {"grid_xt": "x", "grid_yt": "y", "lat": "latitude", "lon": "longitude"}
     )
-    return output
 
-
-def read_gfs_climatology(filename, lat, lon, varname):
-    readin = Dataset(filename)
+    nlev = len(readin.lev.data)
 
     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,
+        DATA = np.empty([nlev, basefile.zc.data.shape[1], basefile.zc.data.shape[2]])
+        for ll in np.arange(nlev):
+            DATA[ll, :, :] = (
+                basefile["zc"].monet.remap_nearest(readin[varname][0, ll, :, :]).data
             )
-
     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 = basefile["zc"].monet.remap_nearest(readin[varname][0, :, :]).data
 
+    readin.close()
     DATA[np.isnan(DATA)] = 0
     DATA[DATA < 0] = 0
     return DATA
 
 
-def read_frp_local(filename, lat, lon, fill_value):
-    readin = Dataset(filename)
+def read_frp_local(filename, basefile):
+    readin = xr.open_dataset(filename)
+    readin = readin.rename({"Latitude": "y", "Longitude": "x"})
+    readin["x"] = readin["x"].where(readin["x"] > 0, readin["x"] + 360)
 
-    # map to met grids
-    xt, yt = np.meshgrid(readin["Longitude"][:], readin["Latitude"][:])
-    xt[xt < 0] = xt[xt < 0] + 360
-    data = np.squeeze(readin["MeanFRP"][:])
+    grid_xt, grid_yt = np.meshgrid(readin["x"].data, readin["y"].data)
+    yt = xr.DataArray(grid_yt, dims=["y", "x"], name="latitude")
+    xt = xr.DataArray(grid_xt, dims=["y", "x"], name="longitude")
+    readin["latitude"] = yt
+    readin["longitude"] = xt
+    readin = readin.set_coords(["latitude", "longitude"])
 
-    DATA = mapping(lat, lon, data.flatten(), yt.flatten(), xt.flatten(), "linear", np.nan)
+    DATA = basefile["zc"].monet.remap_nearest(readin["MeanFRP"][0, :, :]).data
+
+    readin.close()
     return DATA
 
 
@@ -275,24 +256,22 @@ def read_frp_local(filename, lat, lon, fill_value):
     print("------------------------------------")
     print("---- Checking variable dimensions...")
     print("------------------------------------")
-    readin = Dataset(f_met)
-    grid_yt = readin["grid_yt"][:]
-    grid_xt = readin["grid_xt"][:]
-    lat = readin["lat"][:]
-    lon = readin["lon"][:]
-    time = readin["time"][:]
+    basefile = xr.open_dataset(f_met)
+    basefile = basefile.set_coords(["lat", "lon"]).rename(
+        {"grid_xt": "x", "grid_yt": "y", "lat": "latitude", "lon": "longitude"}
+    )
 
     # dimension sizes
-    ntime = len(time)
-    nlat = len(grid_yt)
-    nlon = len(grid_xt)
+    ntime = len(basefile["time"].data)
+    nlat = len(basefile["y"].data)
+    nlon = len(basefile["x"].data)
 
     # var check
-    print("time", time.shape)
-    print("grid_yt", grid_yt.shape)
-    print("grid_xt", grid_xt.shape)
-    print("lat", lat.shape)
-    print("lon", lon.shape)
+    print("time", basefile["time"].data.shape)
+    print("grid_yt", basefile["y"].data.shape)
+    print("grid_xt", basefile["x"].data.shape)
+    print("lat", basefile["latitude"].data.shape)
+    print("lon", basefile["longitude"].data.shape)
 
     """Adding canvar"""
     print("------------------------------------")
@@ -307,27 +286,27 @@ def read_frp_local(filename, lat, lon, fill_value):
         if varname == "lai":
             ATTNAME = ["long_name", "units", "missing_value"]
             ATT = ["Leaf area index", "m^2/m^2", fill_value]
-            DATA = read_gfs_climatology(f_can, lat, lon, "lai")
+            DATA = read_gfs_climatology(f_can, basefile, "lai")
 
         elif varname == "clu":
             ATTNAME = ["long_name", "units", "missing_value"]
             ATT = ["Canopy clumping index", "none", fill_value]
-            DATA = read_gfs_climatology(f_can, lat, lon, "clu")
+            DATA = read_gfs_climatology(f_can, basefile, "clu")
 
         elif varname == "canfrac":
             ATTNAME = ["long_name", "units", "missing_value"]
             ATT = ["Forest fraction of grid cell", "none", fill_value]
-            DATA = read_gfs_climatology(f_can, lat, lon, "canfrac")
+            DATA = read_gfs_climatology(f_can, basefile, "canfrac")
 
         elif varname == "ch":
             ATTNAME = ["long_name", "units", "missing_value"]
             ATT = ["Canopy height above the surface", "m", fill_value]
-            DATA = read_gfs_climatology(f_can, lat, lon, "ch")
+            DATA = read_gfs_climatology(f_can, basefile, "ch")
 
         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")
+            DATA = read_gfs_climatology(f_can, basefile, "pavd")
 
         elif varname == "mol":
             # Reference:
@@ -351,25 +330,28 @@ def read_frp_local(filename, lat, lon, fill_value):
             ATTNAME = ["long_name", "units", "missing_value"]
             ATT = ["Cosine of solar zenith angle", "none", fill_value]
 
+            lat = basefile.latitude.data
+            lon = basefile.longitude.data
+
             time_conv = datetime(
                 int(YY), int(MM), int(DD), int(HH), 0, 0, 0, tzinfo=timezone.utc
             ) + timedelta(hours=int(FH))
             sza = 90 - get_altitude(lat, lon, time_conv)
             DATA = np.cos(sza * 0.0174532925)  # degree to radian
 
-            del [time_conv, sza]
+            del [lat, lon, time_conv, sza]
 
         elif varname == "frp":
             ATTNAME = ["long_name", "units", "missing_value"]
             ATT = ["Mean fire radiative power", "MW", fill_value]
 
             if frp_src == 1:  # 12 month climatology
-                DATA = read_gfs_climatology(f_can, lat, lon, "frp")
+                DATA = read_gfs_climatology(f_can, basefile, "frp")
             elif frp_src == 2:  # ifcanwaf=.FALSE.
                 DATA = np.empty(lat.shape)
                 DATA[:] = 1
             else:
-                DATA = read_frp_local(f_frp, lat, lon, fill_value)
+                DATA = read_frp_local(f_frp, basefile)
 
         elif varname == "href":
             ATTNAME = ["long_name", "units", "missing_value"]