Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add capacity to use AORC forcings via netcdf #16

Draft
wants to merge 13 commits into
base: master
Choose a base branch
from
Draft

Add capacity to use AORC forcings via netcdf #16

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
7b13330
add forcings file creation function to /utils/
Apr 24, 2024
d3b8219
initial commit for nc forcings files
Apr 30, 2024
35a3658
Remove unneeded notebook
May 1, 2024
4641f41
Cleanup
May 1, 2024
633fc62
Refactor aorc function
May 3, 2024
f4d6528
fix import of gdal
May 3, 2024
dd210b7
Fix problem in setting minimum wind forcings
May 3, 2024
d6df23e
Restore ncgs_file
May 3, 2024
bea6f65
Cleanup
May 3, 2024
f6141c2
Don't add to repo. Adding for testing purposes
May 3, 2024
dfd939b
Remove unnecessary import of ncgs_file
May 3, 2024
f24b7ff
Remove unnecessary import
May 3, 2024
2ec8c7f
Remove unnecessary import
May 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
initial commit for nc forcings files
  • Loading branch information
GreyEvenson-NOAA committed Apr 30, 2024
commit d3b8219cc99c75c5f05bedba8323c66f4879b9e9
38 changes: 20 additions & 18 deletions topoflow/components/met_base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@
from topoflow.utils import model_output
from topoflow.utils import rtg_files
from topoflow.utils import time_utils
from topoflow.utils import ncgs_files

#-----------------------------------------------------------------------
class met_component( BMI_base.BMI_component ):
Expand Down Expand Up @@ -1291,10 +1292,10 @@ def update_P_rain(self):
#-------------------------------------------------
P_rain = self.P * (self.T_air > self.T_rain_snow)

if (np.ndim( self.P_rain ) == 0):
if ((np.ndim( self.P_rain ) == 0) & (self.T_air_type.lower() == 'scalar')): # Why is it self.P_rain and not P_rain? LB
self.P_rain.fill( P_rain ) #### (mutable scalar)
else:
self.P_rain[:] = P_rain
self.P_rain = P_rain

if (self.DEBUG):
if (self.P_rain.max() > 0):
Expand Down Expand Up @@ -1344,10 +1345,10 @@ def update_P_snow(self):
#-------------------------------------------------
P_snow = self.P * (self.T_air <= self.T_rain_snow)

if (np.ndim( self.P_snow ) == 0):
if ((np.ndim( self.P_snow ) == 0) & (self.T_air_type.lower() == 'scalar')):
self.P_snow.fill( P_snow ) #### (mutable scalar)
else:
self.P_snow[:] = P_snow
self.P_snow = P_snow

if (self.DEBUG):
if (self.P_snow.max() > 0):
Expand Down Expand Up @@ -2278,8 +2279,9 @@ def read_input_files(self):
else:
units_factor = self.mmph_to_mps
P = model_input.read_next(self.P_unit, self.P_type, rti,
units_factor=units_factor,
NGEN_CSV=self.NGEN_CSV)
units_factor=units_factor,
NGEN_CSV=self.NGEN_CSV,
time_index=self.time_index)

# print('MET: (time,P) =', self.time, P)
# print('shape(P) =', P.shape)
Expand Down Expand Up @@ -2310,7 +2312,7 @@ def read_input_files(self):
# then we'll need to use "fill()" method to prevent breaking
# the reference to the "mutable scalar". (2/7/13)
###############################################################
T_air = model_input.read_next(self.T_air_unit, self.T_air_type, rti)
T_air = model_input.read_next(self.T_air_unit, self.T_air_type, rti, time_index=self.time_index)
## print('## T_air_type =', self.T_air_type )
if (T_air is not None):
# For testing
Expand Down Expand Up @@ -2344,27 +2346,27 @@ def read_input_files(self):
return
###########################################

T_surf = model_input.read_next(self.T_surf_unit, self.T_surf_type, rti)
T_surf = model_input.read_next(self.T_surf_unit, self.T_surf_type, rti, time_index=self.time_index)
if (T_surf is not None):
self.update_var( 'T_surf', T_surf )

RH = model_input.read_next(self.RH_unit, self.RH_type, rti)
RH = model_input.read_next(self.RH_unit, self.RH_type, rti, time_index=self.time_index)
if (RH is not None):
self.update_var( 'RH', RH )

p0 = model_input.read_next(self.p0_unit, self.p0_type, rti)
p0 = model_input.read_next(self.p0_unit, self.p0_type, rti, time_index=self.time_index)
if (p0 is not None):
self.update_var( 'p0', p0 )

uz = model_input.read_next(self.uz_unit, self.uz_type, rti)
uz = model_input.read_next(self.uz_unit, self.uz_type, rti, time_index=self.time_index)
if (uz is not None):
self.update_var( 'uz', uz )

z = model_input.read_next(self.z_unit, self.z_type, rti)
z = model_input.read_next(self.z_unit, self.z_type, rti, time_index=self.time_index)
if (z is not None):
self.update_var( 'z', z )

z0_air = model_input.read_next(self.z0_air_unit, self.z0_air_type, rti)
z0_air = model_input.read_next(self.z0_air_unit, self.z0_air_type, rti, time_index=self.time_index)
if (z0_air is not None):
self.update_var( 'z0_air', z0_air )

Expand All @@ -2374,27 +2376,27 @@ def read_input_files(self):
# Note: We could later write a version of read_next() that takes "self"
# and "var_name" as args and that uses "exec()".
#----------------------------------------------------------------------------
albedo = model_input.read_next(self.albedo_unit, self.albedo_type, rti)
albedo = model_input.read_next(self.albedo_unit, self.albedo_type, rti, time_index=self.time_index)
if (albedo is not None):
self.update_var( 'albedo', albedo )
## self.albedo = albedo

em_surf = model_input.read_next(self.em_surf_unit, self.em_surf_type, rti)
em_surf = model_input.read_next(self.em_surf_unit, self.em_surf_type, rti, time_index=self.time_index)
if (em_surf is not None):
self.update_var( 'em_surf', em_surf )
## self.em_surf = em_surf

dust_atten = model_input.read_next(self.dust_atten_unit, self.dust_atten_type, rti)
dust_atten = model_input.read_next(self.dust_atten_unit, self.dust_atten_type, rti, time_index=self.time_index)
if (dust_atten is not None):
self.update_var( 'dust_atten', dust_atten )
## self.dust_atten = dust_atten

cloud_factor = model_input.read_next(self.cloud_factor_unit, self.cloud_factor_type, rti)
cloud_factor = model_input.read_next(self.cloud_factor_unit, self.cloud_factor_type, rti, time_index=self.time_index)
if (cloud_factor is not None):
self.update_var( 'cloud_factor', cloud_factor )
## self.cloud_factor = cloud_factor

canopy_factor = model_input.read_next(self.canopy_factor_unit, self.canopy_factor_type, rti)
canopy_factor = model_input.read_next(self.canopy_factor_unit, self.canopy_factor_type, rti, time_index=self.time_index)
if (canopy_factor is not None):
self.update_var( 'canopy_factor', canopy_factor )
## self.canopy_factor = canopy_factor
Expand Down
58 changes: 43 additions & 15 deletions topoflow/utils/model_input.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,8 @@
#-------------------------------------------------------------------
import numpy as np
import os.path
import netCDF4 as nc
from . import ncgs_files

#-------------------------------------------------------------------
def open_file(var_type, input_file, NGEN_CSV=False):
Expand Down Expand Up @@ -61,12 +63,21 @@ def open_file(var_type, input_file, NGEN_CSV=False):
#----------------------
line = file_unit.readline()
else:
#--------------------------------------------
# Input file contains a grid or grid stack
# as row-major, binary file with no header.
#--------------------------------------------
file_unit = open(input_file, 'rb')

file_ext = os.path.splitext(input_file)[1]
if file_ext in ('.nc','.nc4'):
#--------------------------------------------
# Input file is netcdf format
#--------------------------------------------
file_unit = nc.Dataset(input_file, mode = 'r')
else:
#--------------------------------------------
# Input file contains a grid or grid stack
# as row-major, binary file with no header.
#--------------------------------------------
file_unit = open(input_file, 'rb')

#print('file_unit:')
#print(file_unit)
return file_unit

# open_file()
Expand Down Expand Up @@ -131,7 +142,7 @@ def read_next2(self, var_name, rti, dtype='float32', factor=1.0):
#-------------------------------------------------------------------
def read_next(file_unit, var_type, rti,
dtype='float32', units_factor=1.0,
NGEN_CSV=False):
NGEN_CSV=False,time_index=None):

if (var_type.lower() == 'scalar'):
#-------------------------------------------
Expand All @@ -157,7 +168,8 @@ def read_next(file_unit, var_type, rti,
# NB! grid_type argument allows DEM to be
# read for GW vars, which might not be FLOAT'
#----------------------------------------------
data = read_grid(file_unit, rti, dtype)
data = read_grid(file_unit, rti, dtype, time_index)

else:
raise RuntimeError('No match found for ' + var_type + '.')
return None
Expand Down Expand Up @@ -223,7 +235,7 @@ def read_scalar(file_unit, dtype='float32',

# read_scalar()
#-------------------------------------------------------------------
def read_grid(file_unit, rti, dtype='float32'):
def read_grid(file_unit, rti, dtype='float32', time_index=None):

#----------------------------------------------------
# Note: Read 2D grid from row-major, binary file.
Expand All @@ -232,13 +244,29 @@ def read_grid(file_unit, rti, dtype='float32'):
# returns an empty array (size 0) with same
# dtype. In this case, return None.
#----------------------------------------------------
grid = np.fromfile(file_unit, count=rti.n_pixels, dtype=dtype)
if (grid.size == 0):
return None

grid = np.reshape(grid, (rti.nrows, rti.ncols))
if (rti.SWAP_ENDIAN):
grid.byteswap(True)
if isinstance(file_unit,nc.Dataset):

#determine name of netcdf variable (from ncgs.get_var_names())
names = list( file_unit.variables.keys() )
dims = list( file_unit.dimensions.keys() )
names = [s for s in names if (s not in dims) and s != 'datetime']
if len(names) != 1:
print('ERROR more than one non dim variables in ',file_unit.filepath())

#get var
grid = file_unit.variables[names[0]][time_index,:,:]

else:

grid = np.fromfile(file_unit, count=rti.n_pixels, dtype=dtype)
if (grid.size == 0):
return None

grid = np.reshape(grid, (rti.nrows, rti.ncols))
if (rti.SWAP_ENDIAN):
grid.byteswap(True)

#-------------------------
return grid

Expand Down
5 changes: 3 additions & 2 deletions topoflow/utils/ncgs_files.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,6 +183,7 @@ def open_file(self, file_name):
#try:
ncgs_unit = nc.Dataset(file_name, mode='r')
self.ncgs_unit = ncgs_unit
self.file_name = file_name
### return ncgs_unit
return True
# except:
Expand Down Expand Up @@ -883,8 +884,8 @@ def get_var_units(self, var_name ):
#----------------------------------------------------------
def get_grid(self, var_name, time_index):

var = self.ncgs_unit.variables[ var_name ]
return var[ time_index ]
var = self.ncgs_unit.variables[var_name][time_index,:,:]
return var

# get_grid()
#----------------------------------------------------------
Expand Down
64 changes: 46 additions & 18 deletions topoflow/utils/regrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@
from osgeo import gdal, osr ## ogr
import glob, sys
import os, os.path
import netCDF4 as nc

from . import rti_files
from . import rtg_files
Expand Down Expand Up @@ -1868,7 +1869,7 @@ def create_rts_from_chirps_files( rts_file='TEST.rts',
# # rename_nc_files()
#-------------------------------------------------------------------

def create_ncgs_forcings_file(var_name=None, nc_file_in=None, ncgs_file_out=None, grid_info=None,
def create_aorc_ncgs_forcings_file(var_name=None, nc_file_in=None, grid_info=None,
resample_algo='bilinear', SILENT=False, VERBOSE=False, NC4=False):

#---------------------------------------------
Expand Down Expand Up @@ -1924,7 +1925,7 @@ def create_ncgs_forcings_file(var_name=None, nc_file_in=None, ncgs_file_out=None
var_units = ncgs_in.get_var_units(var_name=var_name)

#-----------------------------------------
# Create time_info for output nc file from input nc file if none given
# Create time_info for output nc file from input nc file
#-----------------------------------------
time_values = ncgs_in.ncgs_unit.variables['time'][:]
time_dtype = ncgs_in.ncgs_unit.variables['time'].dtype
Expand Down Expand Up @@ -1953,7 +1954,7 @@ def create_ncgs_forcings_file(var_name=None, nc_file_in=None, ncgs_file_out=None
#-----------------------------------------
# Close input nc file
#-----------------------------------------
ncgs_in.close_file() # close file
ncgs_in.close_file()

#-----------------------------------------
# Read var using gdal (to faciliate regridding)
Expand Down Expand Up @@ -2042,39 +2043,67 @@ def create_ncgs_forcings_file(var_name=None, nc_file_in=None, ncgs_file_out=None
print('Converting units: [kg m-2 s-1] to [mmph]...')
w = (grid2 != nc_nodata) # (boolean array)
grid2[w] *= 3600.0 # (preserve nodata)
var_units = 'mmph'
if (var_name == 'Tair_f_inst'):
print('Converting units: Kelvin to Celsius...')
w = (grid2 != nc_nodata) # (boolean array)
grid2[w] -= 273.15 # (preserve nodata)
grid2[w] -= 273.15 # (preserve nodata)
var_units = 'C'
if (var_name == 'Psurf_f_inst'):
print('Converting units: Pa to mbar...')
w = (grid2 != nc_nodata) # (boolean array)
grid2[w] /= 100.0 # (preserve nodata)
var_units = 'mbar'

#-----------------------------------------
# name conversions
#-----------------------------------------
topoflow_var_name = None
if (var_name == 'Rainf_tavg'):
topoflow_var_name = 'P'
elif (var_name == 'Tair_f_inst'):
topoflow_var_name = 'T_air'
elif (var_name == 'Psurf_f_inst'):
topoflow_var_name = 'p0'
elif (var_name == 'Lwnet_tavg'):
topoflow_var_name = 'Qn_LW'
elif (var_name == 'Swnet_tavg'):
topoflow_var_name = 'Qn_SW'
elif (var_name == 'Wind_f_inst'):
topoflow_var_name = 'uz'
if topoflow_var_name is None:
topoflow_var_name = var_name

#-----------------------------------------
# initialize new file
#-----------------------------------------
ncgs_out = ncgs_files.ncgs_file()
ncgs_out.open_new_file(file_name=ncgs_file_out,grid_info=grid_info,time_info=time_info,var_name=var_name,
units_name='m',time_units='minutes',time_res=time_res,
output_file_name = os.path.join(os.path.dirname(nc_file_in),'AORC_'+topoflow_var_name+'.nc')
ncgs_out.open_new_file(file_name=output_file_name,grid_info=grid_info,time_info=time_info,var_name=topoflow_var_name,
units_name=var_units,time_units='minutes',time_res=time_res,
OVERWRITE_OK=True,MAKE_RTI=False)

#------------------------------------
# create attributes for the forcings variable
#------------------------------------
ncgs_out.ncgs_unit.variables[var_name].long_name = var_name
ncgs_out.ncgs_unit.variables[var_name].units = var_units
ncgs_out.ncgs_unit.variables[var_name].dx = grid_info.xres
ncgs_out.ncgs_unit.variables[var_name].dy = grid_info.yres
ncgs_out.ncgs_unit.variables[var_name].y_south_edge = grid_info.y_south_edge
ncgs_out.ncgs_unit.variables[var_name].y_north_edge = grid_info.y_north_edge
ncgs_out.ncgs_unit.variables[var_name].x_west_edge = grid_info.x_west_edge
ncgs_out.ncgs_unit.variables[var_name].x_east_edge = grid_info.x_east_edge
ncgs_out.ncgs_unit.variables[topoflow_var_name].long_name = topoflow_var_name
ncgs_out.ncgs_unit.variables[topoflow_var_name].dx = grid_info.xres
ncgs_out.ncgs_unit.variables[topoflow_var_name].dy = grid_info.yres
ncgs_out.ncgs_unit.variables[topoflow_var_name].y_south_edge = grid_info.y_south_edge
ncgs_out.ncgs_unit.variables[topoflow_var_name].y_north_edge = grid_info.y_north_edge
ncgs_out.ncgs_unit.variables[topoflow_var_name].x_west_edge = grid_info.x_west_edge
ncgs_out.ncgs_unit.variables[topoflow_var_name].x_east_edge = grid_info.x_east_edge

#------------------------------------
# Populate the variable
#------------------------------------
ncgs_out.ncgs_unit.variables[var_name][:,:,:] = grid2
ncgs_out.ncgs_unit.variables[topoflow_var_name][:,:,:] = grid2
ncgs_out.ncgs_unit.variables['time'][:] = [i+1 for i in range(grid2.shape[0])]

#------------------------------------
# Close the file
#------------------------------------
ncgs_out.close()

#---------------------------
# Determine variable units
Expand Down Expand Up @@ -2164,7 +2193,6 @@ def create_ncgs_forcings_file(var_name=None, nc_file_in=None, ncgs_file_out=None
print( 'min(variable) =', vmin, '(possible nodata)' )
print( 'bad_count =', bad_count )
print( 'n_grids =', count )
print( 'Finished saving data to rts file.')
print( ' ')
ncgs_out.close()

# create_ncgs_forcings_file()