Skip to content

Commit

Permalink
erge branch 'main' of github.com:ArcticSnow/TopoPyScale into main
Browse files Browse the repository at this point in the history
  • Loading branch information
@joelfiddes
joelfiddes committed Aug 14, 2024
2 parents 7b35bcd + f5ad251 commit 94e9afb
Show file tree
Hide file tree
Showing 4 changed files with 233 additions and 3 deletions.
167 changes: 165 additions & 2 deletions TopoPyScale/topo_export.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1124,7 +1124,6 @@ def to_crocus(ds,
climate_dataset_name (str): name of original climate dataset. Default 'ERA5',
project_author (str): name of project author(s)
snow_partition_method (str): snow/rain partitioning method: default 'jennings2018_trivariate'
"""
# create one file per point_name
n_digits = len(str(ds.point_name.values.max()))
Expand Down Expand Up @@ -1227,7 +1226,157 @@ def to_crocus(ds,
fo.FORCE_TIME_STEP.attrs = {'units': 's', 'standard_name': 'FORCE_TIME_STEP', 'long_name': 'Forcing_Time_Step',
'_FillValue': -9999999.0}

fo = fo.expand_dims('Number_of_points')
fo = fo.rename_dims({'point_name':'Number_of_points'})
fo = fo.transpose()

fo.to_netcdf(foutput, mode='w', format='NETCDF3_CLASSIC', unlimited_dims={'time': True},
encoding={'time': {'dtype': 'float64', 'calendar': 'standard'},
'Tair': {'dtype': 'float64'},
'Qair': {'dtype': 'float64'},
'Snowf': {'dtype': 'float64'},
'Rainf': {'dtype': 'float64'},
'DIR_SWdown': {'dtype': 'float64'},
'LWdown': {'dtype': 'float64'},
'PSurf': {'dtype': 'float64'},
'HUMREL': {'dtype': 'float64'},
'Wind': {'dtype': 'float64'},
'Wind_DIR': {'dtype': 'float64'},
'CO2air': {'dtype': 'float64'},
'NEB': {'dtype': 'float64'},
'SCA_SWdown': {'dtype': 'float64'},
'ZREF': {'dtype': 'float64'},
'UREF': {'dtype': 'float64'}})
print('---> File {} saved'.format(foutput))


def to_surfex(ds,
df_pts,
fname_format='FORCING*.nc',
scale_precip=1,
climate_dataset_name='ERA5',
project_author='S. Filhol',
snow_partition_method='continuous',
year_start='2023-08-01 00:00:00',
year_end='2024-06-30 23:00:00'):
"""
Function to export toposcale output to Surfex forcing netcdf format. Generates one file per year
Args:
ds (dataset): Toposcale downscaled dataset.
df_pts (dataframe): with point list info (x,y,elevation,slope,aspect,svf,...)
fname_format (str): filename format. point_name is inserted where * is
scale_precip (float): scaling factor to apply on precipitation. Default is 1
climate_dataset_name (str): name of original climate dataset. Default 'ERA5',
project_author (str): name of project author(s)
snow_partition_method (str): snow/rain partitioning method: default 'jennings2018_trivariate'
year_start (str): start of year for yearly Surfex forcing
ToDo:
- [x] split to yearly files,
"""
# create one file per point_name
n_digits = len(str(ds.point_name.max().values))
ver_dict = tu.get_versionning()

tvec_start, tvec_end = tu.time_vecs(start_date=year_start, end_date=year_end)

for year_id, start in enumerate(tvec_start):
fo = ds.sel(time=slice(start, tvec_end[year_id])).copy()
start_str = start.strftime('%Y%m%d%H')
end_str = tvec_end[year_id].strftime('%Y%m%d%H')
foutput = f"{str(fname_format).split('*')[0]}_{start_str}06_{end_str}06.nc"

fo['point_name'] = fo.point_name.astype(int)
fo = fo.rename_dims({'point_name':'Number_of_points'})
fo = fo.rename({'t':'Tair',
'SW':'DIR_SWdown',
'LW':'LWdown',
'p':'PSurf',
'q':'Qair',
'ws':'Wind'})
fo['Wind_DIR'] = fo.wd * 180 / np.pi
fo['HUMREL'] = (('Number_of_points','time'), mu.q_2_rh(fo.Tair.values, fo.PSurf.values, fo.Qair.values) * 100)
fo['precip'] = (('Number_of_points','time'), fo.tp.values / 3600 * scale_precip) # convert from mm/hr to mm/s
rh = mu.q_2_rh(fo.Tair.values, fo.PSurf.values, fo.Qair.values)
Rainf, Snowf = mu.partition_snow(fo.precip.values, fo.Tair.values, rh, fo.PSurf.values,
method=snow_partition_method)
fo['Rainf'] = (('Number_of_points','time'), Rainf)
fo['Snowf'] = (('Number_of_points','time'), Snowf)

fo['NEB'] = (('Number_of_points','time'), np.zeros(fo.Tair.shape))
fo['CO2air'] = (('Number_of_points','time'), np.zeros(fo.Tair.shape))
fo['SCA_SWdown'] = (('Number_of_points','time'), np.zeros(fo.Tair.shape))

fo = fo.drop_vars(['precip', 'u', 'v', 'w', 'wd', 'vp', 'cse','SW_direct', 'cos_illumination', 'SW_diffuse', 'precip_lapse_rate', 'tp' ])

fo.Tair.attrs = {'units': 'K', 'standard_name': 'Tair', 'long_name': 'Near Surface Air Temperature',
'_FillValue': -9999999.0}
fo.Qair.attrs = {'units': 'kg/kg', 'standard_name': 'Qair', 'long_name': 'Near Surface Specific Humidity',
'_FillValue': -9999999.0}
fo.HUMREL.attrs = {'units': '%', 'standard_name': 'HUMREL', 'long_name': 'Relative Humidity',
'_FillValue': -9999999.0}
fo.Wind.attrs = {'units': 'm/s', 'standard_name': 'Wind', 'long_name': 'Wind Speed', '_FillValue': -9999999.0}
fo.Wind_DIR.attrs = {'units': 'deg', 'standard_name': 'Wind_DIR', 'long_name': 'Wind Direction',
'_FillValue': -9999999.0}
fo.Rainf.attrs = {'units': 'kg/m2/s', 'standard_name': 'Rainf', 'long_name': 'Rainfall Rate',
'_FillValue': -9999999.0}
fo.Snowf.attrs = {'units': 'kg/m2/s', 'standard_name': 'Snowf', 'long_name': 'Snowfall Rate',
'_FillValue': -9999999.0}
fo.DIR_SWdown.attrs = {'units': 'W/m2', 'standard_name': 'DIR_SWdown',
'long_name': 'Surface Incident Direct Shortwave Radiation', '_FillValue': -9999999.0}
fo.LWdown.attrs = {'units': 'W/m2', 'standard_name': 'LWdown',
'long_name': 'Surface Incident Longtwave Radiation', '_FillValue': -9999999.0}
fo.PSurf.attrs = {'units': 'Pa', 'standard_name': 'PSurf', 'long_name': 'Surface Pressure',
'_FillValue': -9999999.0}
fo.SCA_SWdown.attrs = {'units': 'W/m2', 'standard_name': 'SCA_SWdown',
'long_name': 'Surface Incident Diffuse Shortwave Radiation', '_FillValue': -9999999.0}
fo.CO2air.attrs = {'units': 'kg/m3', 'standard_name': 'CO2air', 'long_name': 'Near Surface CO2 Concentration',
'_FillValue': -9999999.0}
fo.NEB.attrs = {'units': 'between 0 and 1', 'standard_name': 'NEB', 'long_name': 'Nebulosity',
'_FillValue': -9999999.0}

fo.attrs = {'title': 'Forcing for SURFEX CROCUS',
'source': 'data from {} downscaled with TopoPyScale ready for CROCUS'.format(climate_dataset_name),
'creator_name': 'Dataset created by {}'.format(project_author),
'package_version': ver_dict.get('package_version'),
'git_commit': ver_dict.get('git_commit'),
'url_TopoPyScale': 'https://github.com/ArcticSnow/TopoPyScale',
'date_created': dt.datetime.now().strftime('%Y/%m/%d %H:%M:%S')}

ds_pts = xr.Dataset(df_pts).rename_dims({'dim_0':'Number_of_points'})
fo['ZS'] = (('Number_of_points'), np.float64(ds_pts.elevation.sel(Number_of_points=fo.Number_of_points).values))
fo['aspect'] = (('Number_of_points'), np.float64(np.rad2deg(ds_pts.aspect.sel(Number_of_points=fo.Number_of_points).values)))
fo['slope'] = (('Number_of_points'), np.float64(np.rad2deg(ds_pts.slope.sel(Number_of_points=fo.Number_of_points).values)))
fo['massif_number'] = (('Number_of_points'), np.zeros(fo.ZS.shape) + np.float64(0))
fo['LAT'] = (('Number_of_points'), ds_pts.latitude.sel(Number_of_points=fo.Number_of_points).values)
fo['LON'] = (('Number_of_points'), ds_pts.longitude.sel(Number_of_points=fo.Number_of_points).values)
ds_pts=None
fo.LAT.attrs = {'units': 'degrees_north', 'standard_name': 'LAT', 'long_name': 'latitude',
'_FillValue': -9999999.0}
fo.LON.attrs = {'units': 'degrees_east', 'standard_name': 'LON', 'long_name': 'longitude',
'_FillValue': -9999999.0}
fo['ZREF'] = (('Number_of_points'), np.zeros(fo.ZS.shape) + np.float64(2))
fo['UREF'] = (('Number_of_points'), np.zeros(fo.ZS.shape) + np.float64(10))
fo.ZREF.attrs = {'units': 'm', 'standard_name': 'ZREF', 'long_name': 'Reference_Height',
'_FillValue': -9999999.0}
fo.UREF.attrs = {'units': 'm', 'standard_name': 'UREF', 'long_name': 'Reference_Height_for_Wind',
'_FillValue': -9999999.0}
fo.aspect.attrs = {'units': 'degree from north', 'standard_name': 'aspect', 'long_name': 'slope aspect',
'_FillValue': -9999999.0}
fo.slope.attrs = {'units': 'degrees from horizontal', 'standard_name': 'slope', 'long_name': 'slope angle',
'_FillValue': -9999999.0}
fo.ZS.attrs = {'units': 'm', 'standard_name': 'ZS', 'long_name': 'altitude', '_FillValue': -9999999.0}
fo.massif_number.attrs = {'units': '', 'standard_name': 'massif_number',
'long_name': 'SAFRAN massif number in SOPRANO world', '_FillValue': -9999999.0}
fo.time.attrs = {'standard_name': 'time', 'long_name': 'time', '_FillValue': -9999999.0}
fo['FRC_TIME_STP'] = (('Number_of_points'), np.zeros(fo.ZS.shape) + np.float64(3600))
fo.FRC_TIME_STP.attrs = {'units': 's', 'standard_name': 'FRC_TIME_STP', 'long_name': 'Forcing_Time_Step',
'_FillValue': -9999999.0}
#fo['FORCE_TIME_STEP'] = (('Number_of_points'), np.zeros(fo.ZS.shape) + np.float64(3600))
#fo.FORCE_TIME_STEP.attrs = {'units': 's', 'standard_name': 'FORCE_TIME_STEP', 'long_name': 'Forcing_Time_Step',
# '_FillValue': -9999999.0}

fo = fo.transpose()

fo.to_netcdf(foutput, mode='w', format='NETCDF3_CLASSIC', unlimited_dims={'time': True},
Expand All @@ -1248,6 +1397,20 @@ def to_crocus(ds,
'ZREF': {'dtype': 'float64'},
'UREF': {'dtype': 'float64'}})
print('---> File {} saved'.format(foutput))
# clear memory
fo = None














def to_snowpack(ds, fname_format='smet_pt_*.tx'):
Expand Down
40 changes: 40 additions & 0 deletions TopoPyScale/topo_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,5 +300,45 @@ def plot_xyline(ax):
#FsmPlot(df)


def time_vecs(start_date, end_date):
"""Function to build timestamp vectors of year starts and year ends for arbitrary time range.
Args:
start_date (str): _description_
end_date (str): _description_
Returns:
start_year: vector of start of years (datetime64)
end_year: vector of end of year (datetime64)
"""

current_start = pd.to_datetime(start_date)
end_date = pd.to_datetime(end_date)
start_year = []
end_year = []

while current_start <= end_date:
# Append the current start date
start_year.append(current_start)

# Calculate the end date for the current cycle
next_start = current_start + pd.DateOffset(years=1)
current_end = next_start - pd.DateOffset(days=1)

# If the cycle's end date goes beyond the overall end date, break and add the final segment
if current_end > end_date:
current_end = end_date
end_year.append(current_end)
break

# Append the current end date
end_year.append(current_end)

# Update current start date to next year start
current_start = next_start

start_year = pd.to_datetime(start_year)
end_year = pd.to_datetime(end_year)

return start_year, end_year

25 changes: 25 additions & 0 deletions TopoPyScale/topoclass.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -780,6 +780,31 @@ def to_crocus(self, fname_format='CROCUS_pt_*.nc', scale_precip=1):
climate_dataset_name=self.config.project.climate,
project_author=self.config.project.authors)

def to_surfex(self,
fname_format='FORCING_*.nc',
scale_precip=1,
year_start=None,
year_end=None):
"""
function to export toposcale output to surfex format .nc. This functions saves one file per point_name
Args:
fout_format (str): filename format. point_name is inserted where * is
scale_precip(float): scaling factor to apply on precipitation. Default is 1
"""
if year_start is None or year_end is None:
raise ValueError("Arguments year_start and year_end are required.")

te.to_surfex(self.downscaled_pts,
self.toposub.df_centroids,
fname_format=self.config.outputs.path / fname_format,
scale_precip=scale_precip,
climate_dataset_name=self.config.project.climate,
project_author=self.config.project.authors,
snow_partition_method='continuous',
year_start=year_start,
year_end=year_end)

def to_snowmodel(self, fname_format='Snowmodel_stn_*.csv'):
"""
function to export toposcale output to snowmodel format .ascii, for single station standard
Expand Down
4 changes: 3 additions & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@
packages=find_packages(),
python_requires='>=3.8',
install_requires=['xarray[complete]',
'rioxarray',
'matplotlib',
'scikit-learn',
'pandas',
Expand All @@ -69,6 +70,7 @@
'pyproj',
'dask',
'configobj',
'munch'],
'munch',
'gitpython'],
include_package_data=True
)

0 comments on commit 94e9afb

Please sign in to comment.