Skip to content

Commit

Permalink
Update global_data_process.py
Browse files Browse the repository at this point in the history
  • Loading branch information
@angehung5
angehung5 authored Nov 7, 2023
1 parent 6891404 commit 8ecd4a7
Showing 1 changed file with 23 additions and 51 deletions.
74 changes: 23 additions & 51 deletions python/global_data_process.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,15 @@
from scipy.interpolate import griddata

"""User Arguments"""
## Time: yyyymmddhhfff
# Time: yyyymmddhhfff
timelist = sys.argv[1]
timelist = np.array(timelist.split(",")).astype(str)


"""User Options"""
path = "/scratch/pcampbe8/canopy-app/input" # work directory
ref_lev = 10 # reference height (m, a.g.l.)
frp_src = 0 # frp data source (0: local source; 1: check local source first, switch to climatological file if no available data; 2: 12 month climatology; 3: all ones when ifcanwaf=.FALSE.)
path = "/scratch/pcampbe8/canopy-app/input" # work directory
ref_lev = 10 # reference height (m, a.g.l.)
frp_src = 0 # frp data source (0: local source; 1: check local source first, switch to climatological file if no available data; 2: 12 month climatology; 3: all ones when ifcanwaf=.FALSE.)



Expand All @@ -47,7 +47,6 @@
# ------------------------------------------------------------------------- #



"""Global Settings"""
# domain
lat_lim = [-90, 90]
Expand Down Expand Up @@ -89,30 +88,33 @@ def read_varatt(var):
att = [var.getncattr(X) for X in attname]
return attname, att


def write_varatt(var, attname, att):
for X in np.arange(len(attname)):
if attname[X] == "_FillValue":
continue
elif attname[X] == "missing_value":
value = np.float32(att[X])
value = np.round(value / 1e+15) * 1e+15
value = np.round(value / 1e15) * 1e15
var.setncattr(attname[X], value)
else:
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
)
return output


def read_gfs_climatology(filename, lat, lon, varname):
readin = Dataset(filename)

if varname == "pavd":
# map to met grids
yt = readin["lat"][:]
xt = readin["lon"][:]
yt = readin["lat"][:]
xt = readin["lon"][:]
data = np.squeeze(readin[varname][:])

DATA = np.empty([data.shape[0], lat.shape[0], lat.shape[1]])
Expand All @@ -125,13 +127,13 @@ def read_gfs_climatology(filename, lat, lon, varname):
yt.flatten(),
xt.flatten(),
"linear",
np.nan
np.nan,
)

else:
# map to met grids
yt = readin["lat"][:]
xt = readin["lon"][:]
yt = readin["lat"][:]
xt = readin["lon"][:]
data = np.squeeze(readin[varname][0, :, :])

DATA = mapping(
Expand All @@ -142,12 +144,13 @@ def read_gfs_climatology(filename, lat, lon, varname):
DATA[DATA < 0] = 0
return DATA


def read_frp_local(filename, lat, lon, fill_value):
readin = Dataset(filename)

# map to met grids
xt, yt = np.meshgrid(readin["Longitude"][:], readin["Latitude"][:])
xt[xt < 0] = xt[xt < 0]+360
xt[xt < 0] = xt[xt < 0] + 360
data = np.squeeze(readin["MeanFRP"][:])

DATA = mapping(
Expand All @@ -156,18 +159,13 @@ def read_frp_local(filename, lat, lon, fill_value):
return DATA




starttime = datetime.now()
print("------------------------------------")
print("---- Global input pre-process start!", starttime.strftime("%Y/%m/%d %H:%M:%S"))
print("------------------------------------")




for inputtime in timelist:

if (
(len(inputtime) != 13)
or (int(inputtime[4:6]) > 12)
Expand All @@ -181,7 +179,6 @@ def read_frp_local(filename, lat, lon, fill_value):
print("---- " + inputtime[:10] + " f" + inputtime[10:] + " processing...")
print("------------------------------------")


"""Settings"""
# date/time
YY = inputtime[:4] # year
Expand All @@ -190,7 +187,6 @@ def read_frp_local(filename, lat, lon, fill_value):
HH = inputtime[8:10] # hour in UTC
FH = inputtime[10:] # forecast hour, for met file only


# input/output files
f_met = (
path + "/gfs.t" + HH + "z." + YY + MM + DD + ".sfcf" + FH + ".nc"
Expand Down Expand Up @@ -222,8 +218,6 @@ def read_frp_local(filename, lat, lon, fill_value):
elif frp_src == 2: # climatological frp
f_frp = path + "/gfs.canopy.t" + HH + "z." + YY + MM + DD + ".sfcf000.nc"



"""Data Download"""
"""Download from servers if required files do not exist."""
print("---- Checking required files...")
Expand Down Expand Up @@ -313,7 +307,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" + HH + "z." + YY + MM + DD + ".sfcf000.nc"
f_frp = path + "/gfs.canopy.t" + HH + "z." + YY + MM + DD + ".sfcf000.nc"

if frp_src == 2: # 12 month climatology frp
if os.path.isfile(f_frp) == True:
Expand Down Expand Up @@ -346,41 +340,31 @@ def read_frp_local(filename, lat, lon, fill_value):
print("-----------!!!WARNING!!!-------------")
print("---!!!Climatological FRP is used!!!--")




os.system("cp " + f_met + " " + f_output) # copy gfs met file for appending




"""Reading dimensions"""
print("------------------------------------")
print("---- Checking variable dimensions...")
print("------------------------------------")
readin = Dataset(f_met)
readin = Dataset(f_met)
grid_yt = readin["grid_yt"][:]
grid_xt = readin["grid_xt"][:]
lat = readin["lat"][:]
lon = readin["lon"][:]
time = readin["time"][:]


# dimension sizes
ntime = len(time)
nlat = len(grid_yt)
nlon = len(grid_xt)


# 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)



"""Adding canvar"""
print("------------------------------------")
print("---- Generating canopy variables...")
Expand All @@ -396,31 +380,26 @@ def read_frp_local(filename, lat, lon, fill_value):
ATT = ["Leaf area index", "m^2/m^2", fill_value]
DATA = read_gfs_climatology(f_can, lat, lon, "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")


elif varname == "ffrac":
ATTNAME = ["long_name", "units", "missing_value"]
ATT = ["Forest fraction of grid cell", "none", fill_value]
DATA = read_gfs_climatology(f_can, lat, lon, "ffrac")


elif varname == "fh":
ATTNAME = ["long_name", "units", "missing_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
Expand All @@ -434,12 +413,11 @@ def read_frp_local(filename, lat, lon, fill_value):
shtfl = np.squeeze(readin["shtfl"][:])

DATA = (-1 * den * Cp * t2m * (fricv**3)) / (K * g * shtfl)
DATA[DATA > 500] = 500
DATA[DATA > 500] = 500
DATA[DATA <- 500] = -500

del [readin, t2m, fricv, shtfl]


elif varname == "csz":
ATTNAME = ["long_name", "units", "missing_value"]
ATT = ["Cosine of solar zenith angle", "none", fill_value]
Expand All @@ -451,8 +429,7 @@ def read_frp_local(filename, lat, lon, fill_value):
DATA = np.cos(sza * 0.0174532925) # degree to radian

del [time_conv, sza]



elif varname == "frp":
ATTNAME = ["long_name", "units", "missing_value"]
ATT = ["Mean fire radiative power", "MW", fill_value]
Expand All @@ -465,29 +442,26 @@ def read_frp_local(filename, lat, lon, fill_value):
else:
DATA = read_frp_local(f_frp, lat, lon, fill_value)


elif varname == "href":
ATTNAME = ["long_name", "units", "missing_value"]
ATT = ["Reference height above the surface", "m", fill_value]
DATA = np.empty([nlat, nlon])
DATA[:] = ref_lev


# var check
print("Dimension/Attributes:")
print(DATA.shape)
print(ATTNAME)
print(ATT)


# adding to output file
if varname == "pavd":
output = Dataset(f_output, "a")
output.createDimension("level", DATA.shape[0])

var_lev = output.createVariable("lev", "float", ("level", ))
var_bot = output.createVariable("layer_bottom", "i4", ("level", ))
var_top = output.createVariable("layer_top", "i4", ("level", ))
var_lev = output.createVariable("lev", "float", ("level",))
var_bot = output.createVariable("layer_bottom", "i4", ("level",))
var_top = output.createVariable("layer_top", "i4", ("level",))
var = output.createVariable(
varname,
"float",
Expand Down Expand Up @@ -541,11 +515,9 @@ def read_frp_local(filename, lat, lon, fill_value):
print("------------------------------------")



endtime = datetime.now()



print("---- Global input pre-process complete!", endtime.strftime("%Y/%m/%d %H:%M:%S"))
print("---- Process time:", str(endtime-starttime))
print("---- Process time:", str(endtime - starttime))
print("------------------------------------")

0 comments on commit 8ecd4a7

Please sign in to comment.