diff --git a/cosipy/utilities/aws2cosipy/aws2cosipy.py b/cosipy/utilities/aws2cosipy/aws2cosipy.py index 2e541e9c..1b1a83c5 100644 --- a/cosipy/utilities/aws2cosipy/aws2cosipy.py +++ b/cosipy/utilities/aws2cosipy/aws2cosipy.py @@ -35,6 +35,7 @@ --xr Right longitude value of the subset. --yl Lower latitude value of the subset. --yu Upper latitude value of the subset. + --sw Path to the HORAYZON LUT. """ @@ -451,6 +452,7 @@ def create_2D_input( x1=None, y0=None, y1=None, + corr_file=None ): """Create a 2D input dataset from a .csv file. @@ -646,6 +648,41 @@ def create_2D_input( ) else: G_interp[t, i, j] = sw[t] + + elif _cfg.radiation["radiationModule"] == "Horayzon2022": + print("Run the radiation moduel HORAYZON") + + #get correction factor which must be computed beforehand + try: + correction_factor = xr.open_dataset(corr_file) + except: + print("There is no HORAYZON LUT loaded. Please ensure you parsed the correct path.") + raise_nan_error() + + #impose limits on correction factor to ensure reasonable range + corr_vals = correction_factor['sw_dir_cor'].values + corr_vals_clip = np.where(corr_vals > 25, 25, corr_vals) + correction_factor['sw_dir_cor'] = (('time', 'lat', 'lon'), corr_vals_clip) + correction_factor['doy'] = correction_factor.time.dt.dayofyear + correction_factor['hour'] = correction_factor.time.dt.hour + #Create unique identifier for doy and hour and set as index + correction_factor['time_id'] = correction_factor['doy'] + correction_factor['hour']/100 + correction_factor = correction_factor.set_index(time="time_id") + + #Start loop over each timestep + for t in range(time_index): + #get doy and timestep from df and check with doy and hour from static file + year = df.index[t].year + doy = df.index[t].dayofyear + hour = df.index[t].hour + #we rely on the fact that HORAYZON LUT was created for a leap year and factors dont change with time + if (year % 4 != 0 and doy > 59): #59th DOY = February 28th, leap year continues with 29th + doy = doy + 1 + + time_identifier = doy + hour/100 + sw_cor_val = correction_factor.sel(time=time_identifier)['sw_dir_cor'] + #multiply correction factors with forcing + G_interp[t,:,:] = sw_cor_val * sw[t] elif _cfg.radiation["radiationModule"] == "Moelg2009": print("Run the radiation module Moelg2009") @@ -1001,6 +1038,14 @@ def get_user_arguments(parser: argparse.ArgumentParser) -> argparse.Namespace: const=None, help="Upper latitude value of the subset", ) + parser.add_argument( + "--sw", + dest="corr_file", + type=str, + metavar="", + const=None, + help="Path to the HORAYZON LUT table", + ) arguments = parser.parse_args() return arguments @@ -1047,6 +1092,7 @@ def main(): _args.xr, _args.yl, _args.yu, + _args.corr_file ) diff --git a/cosipy/utilities/aws2cosipy/createHORAYZONLUT.py b/cosipy/utilities/aws2cosipy/createHORAYZONLUT.py new file mode 100644 index 00000000..6e5ff940 --- /dev/null +++ b/cosipy/utilities/aws2cosipy/createHORAYZONLUT.py @@ -0,0 +1,672 @@ +""" +Creates the ray-tracing based correction factors for incoming shortwave radiation. +Processing routine is based on the HORAYZON package (https://github.com/ChristianSteger/HORAYZON). +Please refer to the publication by Steger et al., (2022): https://gmd.copernicus.org/articles/15/6817/2022/. + +This script requires the installation of HORAYZON and xesmf (https://xesmf.readthedocs.io/en/stable/). + +Prerequisite is a static file that incorporates surrounding terrain, created through the COSIPY utilities. +If regridding is desired, correction factors should be calculated at high resolution first and then regridded to lower resolution +instead of calculation correction factors on low resolution static files directly. + +Usage: + +python -m cosipy.utilities.aws2cosipy.createHORAYZONLUT -s -o [-c ]\ + [-r ] [-e ] [-es ] [-d Path to the static file that should be used for the calculation. + -o, --output Path to the resulting netCDF file. + +All other arguments are supplied with a default behaviour of -r False, -e False, -c None, -es 30, -d None +Elevation bands are only calculated when regridding is inactive which will be enforced if -e True. +Optional arguments: + -c, --coarse-static Path to the coarse static file to which the results should be regridded. Only used when -r True + -r, --regridding String conveted to boolean information on whether to regrid or not + -e, --elevation_prof String converted to boolean information on whether to compute 1D elevation bands or not + -es, --elevation_size Elevation bin size, used when -e True + -d, --elev_data Only used when -e True: Path where elevation bands static file should be written to. +""" + +import argparse +import numpy as np +from skyfield.api import load, wgs84 +import time +import datetime as dt +import horayzon as hray +import xarray as xr +import xesmf as xe +import pandas as pd + +from cosipy.utilities.config_utils import UtilitiesConfig + +_args = None +_cfg = None + +ellps = "WGS84" # Earth's surface approximation (sphere, GRS80 or WGS84) + +# ---------------------------- +# Some helper functions +# ---------------------------- + +#https://stackoverflow.com/a/43357954 +def str2bool(v): + if isinstance(v, bool): + return v + if v.lower() in ('yes', 'true', 't', 'y', '1'): + return True + elif v.lower() in ('no', 'false', 'f', 'n', '0'): + return False + else: + raise argparse.ArgumentTypeError('Boolean value expected.') + +def add_variable_along_timelatlon(ds, var, name, units, long_name): + """ This function adds missing variables to the DATA class """ + ds[name] = (('time','lat','lon'), var) + ds[name].attrs['units'] = units + ds[name].attrs['long_name'] = long_name + return ds + +def add_variable_along_latlon(ds, var, name, units, long_name): + """ This function self.adds missing variables to the self.DATA class """ + ds[name] = (('lat','lon'), var) + ds[name].attrs['units'] = units + ds[name].attrs['long_name'] = long_name + ds[name].encoding['_FillValue'] = -9999 + return ds + +### function to assign attributes to variable +def assign_attrs(ds, name, units, long_name): + ds[name].attrs['units'] = units + ds[name].attrs['long_name'] = long_name + ds[name].attrs['_FillValue'] = -9999 + +#function for mean of circular values strongly inspired by http://webspace.ship.edu/pgmarr/Geo441/Lectures/Lec%201 +def aspect_means(x): + mean_sine = np.nanmean(np.sin(np.radians(x))) + mean_cosine = np.nanmean(np.cos(np.radians(x))) + r = np.sqrt(mean_cosine**2 + mean_sine**2) + cos_mean = mean_cosine/r + sin_mean = mean_sine/r + mean_angle = np.arctan2(sin_mean, cos_mean) + return np.degrees(mean_angle) + + +def calculate_1d_elevationband(xds, elevation_var, mask_var, var_of_interest, elev_bandsize, slice_idx=None): + + ## first mask vals + xds = xds.where(xds[mask_var] == 1, drop=True) + + #test groupby bins + full_elev_range = xds[elevation_var].values[xds[mask_var] == 1] + bins = np.arange(np.nanmin(full_elev_range), np.nanmax(full_elev_range)+elev_bandsize, elev_bandsize) + labels = bins[:-1] + elev_bandsize/2 + + if var_of_interest in ["lat","lon"]: + values = [] + for i in (bins): + sub = xds.where((xds[mask_var] == 1) & (xds[elevation_var] >= i) & (xds[elevation_var] < i + elev_bandsize), drop=True) + values.append(np.nanmean(sub[var_of_interest].values)) + elif var_of_interest == "ASPECT": + elvs = xds[elevation_var].values.flatten()[xds[mask_var].values.flatten() == 1] + aspects = xds[var_of_interest].values.flatten()[xds[mask_var].values.flatten() == 1] + values = [] + for i in (bins): + values.append(aspect_means(aspects[np.logical_and(elvs >= i, elvs < i+elev_bandsize)])) + elif var_of_interest == mask_var: + if slice_idx is None: + values = xds[var_of_interest].groupby_bins(xds[elevation_var], bins, labels=labels, include_lowest=True).sum(skipna=True, min_count=1) + else: + values = xds[var_of_interest][slice_idx].groupby_bins(xds[elevation_var][slice_idx], bins, labels=labels, include_lowest=True).sum(skipna=True, min_count=1) + ## below calculation doesnt work + #elif var_of_interest in ["aspect","ASPECT"]: + # if slice_idx is None: + # values = xds[var_of_interest].groupby_bins(xds[elevation_var], bins, labels=labels, include_lowest=True).map(aspect_means) + # else: + # values = xds[var_of_interest][slice_idx].groupby_bins(xds[elevation_var][slice_idx], bins, labels=labels, include_lowest=True).map(aspect_means) + else: + if slice_idx is None: + values = xds[var_of_interest].groupby_bins(xds[elevation_var], bins, labels=labels, include_lowest=True).mean(skipna=True) + + else: + values = xds[var_of_interest][slice_idx].groupby_bins(xds[elevation_var][slice_idx], bins, labels=labels, include_lowest=True).mean(skipna=True) + + return values + +def construct_1d_dataset(df): + elev_ds = df.to_xarray() + elev_ds.lon.attrs['standard_name'] = 'lon' + elev_ds.lon.attrs['long_name'] = 'longitude' + elev_ds.lon.attrs['units'] = 'Average Lon of elevation bands' + + elev_ds.lat.attrs['standard_name'] = 'lat' + elev_ds.lat.attrs['long_name'] = 'latitude' + elev_ds.lat.attrs['units'] = 'Average Lat of elevation bands' + assign_attrs(elev_ds, 'HGT','meters','Mean of elevation range per bin as meter above sea level') + assign_attrs(elev_ds, 'ASPECT','degrees','Mean Aspect of slope') + assign_attrs(elev_ds, 'SLOPE','degrees','Mean Terrain slope') + assign_attrs(elev_ds, 'MASK','boolean','Glacier mask') + assign_attrs(elev_ds, 'N_Points','count','Number of Points in each bin') + assign_attrs(elev_ds, 'sw_dir_cor','-','Average shortwave radiation correction factor per elevation band') + + return elev_ds + +def compute_and_slice(latitudes, longitudes, mask_obj): + # Compute indices of inner domain -> needs to encompass everything in range for aggregation + slice_in = (slice(1,latitudes.shape[0]-1, None), slice(1, longitudes.shape[0]-1)) + + # Compute glacier mask + mask_glacier_original = mask_obj + #set NaNs to zero, relict from create static file + mask_glacier_original[np.isnan(mask_glacier_original)] = 0 + mask_glacier = mask_glacier_original.astype(bool) + mask_glacier = mask_glacier[slice_in] #-1 -1 verywhere + + #mask with buffer for aggregation to lower spatial resolutions + #set +- 11 grid cells to "glacier" to allow ensure regridding + ilist = [] + jlist = [] + ## Note that this list is not based on the original shape, see slice_in above + for i in np.arange(0,mask_glacier.shape[0]): + for j in np.arange(0,mask_glacier.shape[1]): + if mask_glacier[i,j] == True: + #print("Grid cell is glacier.") + ilist.append(i) + jlist.append(j) + #create buffer around glacier + ix_latmin = np.min(ilist) + ix_latmax = np.max(ilist) + ix_lonmin = np.min(jlist) + ix_lonmax = np.max(jlist) + + #Watch out that the large domain incorporates the buffer - here selected 11 grid cells + slice_buffer = (slice(ix_latmin-11,ix_latmax+11), slice(ix_lonmin-11, ix_lonmax+11)) + mask_glacier[slice_buffer] = True + + return (slice_in, slice_buffer, mask_glacier, mask_glacier_original) + +def compute_coords(lat, lon, elevation, slice_in): + # Compute ECEF coordinates + x_ecef, y_ecef, z_ecef = hray.transform.lonlat2ecef(*np.meshgrid(lon, lat), + elevation, ellps=ellps) + + # Compute ENU coordinates + trans_ecef2enu = hray.transform.TransformerEcef2enu( + lon_or=lon[int(len(lon) / 2)], lat_or=lat[int(len(lat) / 2)], ellps=ellps) + x_enu, y_enu, z_enu = hray.transform.ecef2enu(x_ecef, y_ecef, z_ecef, + trans_ecef2enu) + + # Compute unit vectors (up and north) in ENU coordinates for inner domain + vec_norm_ecef = hray.direction.surf_norm(*np.meshgrid(lon[slice_in[1]], + lat[slice_in[0]])) + vec_north_ecef = hray.direction.north_dir(x_ecef[slice_in], y_ecef[slice_in], + z_ecef[slice_in], vec_norm_ecef, + ellps=ellps) + del x_ecef, y_ecef, z_ecef + vec_norm_enu = hray.transform.ecef2enu_vector(vec_norm_ecef, trans_ecef2enu) + vec_north_enu = hray.transform.ecef2enu_vector(vec_north_ecef, trans_ecef2enu) + del vec_norm_ecef, vec_north_ecef + + # Merge vertex coordinates and pad geometry buffer + # holds all the data + vert_grid = hray.auxiliary.rearrange_pad_buffer(x_enu, y_enu, z_enu) + + # Compute rotation matrix (global ENU -> local ENU) + + rot_mat_glob2loc = hray.transform.rotation_matrix_glob2loc(vec_north_enu, + vec_norm_enu) + + del vec_north_enu + + # Compute slope (in global ENU coordinates!) + slice_in_a1 = (slice(slice_in[0].start - 1, slice_in[0].stop + 1), + slice(slice_in[1].start - 1, slice_in[1].stop + 1)) + + ## Slope vs plain method -> for comparison later + vec_tilt_enu = \ + np.ascontiguousarray(hray.topo_param.slope_vector_meth( + x_enu[slice_in_a1], y_enu[slice_in_a1], z_enu[slice_in_a1], + rot_mat=rot_mat_glob2loc, output_rot=False)[1:-1, 1:-1]) + + return (vert_grid, vec_norm_enu, vec_tilt_enu, trans_ecef2enu, + x_enu, y_enu, z_enu, rot_mat_glob2loc) + +def merge_timestep_files(datasets, regrid, ds_coarse, static_ds, + elevation_profile, mask_glacier_original, slice_in, slice_buffer, + ): + #Merge single timestep files + now = time.time() + ds_sw_cor = xr.concat(datasets, dim='time') + ds_sw_cor['time'] = pd.to_datetime(ds_sw_cor['time'].values) + if regrid == True: + ds_sw_cor['MASK'] = ds_coarse["MASK"] #replace with original mask, should have same dimensions + else: + #dont have same dimensions + if elevation_profile == True: + ds_sw_cor['HGT'] = ds_sw_cor['HGT'].isel(time=0) + ds_sw_cor['ASPECT'] = ds_sw_cor['ASPECT'].isel(time=0) + ds_sw_cor['SLOPE'] = ds_sw_cor['SLOPE'].isel(time=0) + ds_sw_cor['MASK'] = ds_sw_cor['MASK'].isel(time=0) + ds_sw_cor['N_Points'] = ds_sw_cor['N_Points'].isel(time=0) + else: + mask_holder = mask_glacier_original[slice_in] + add_variable_along_latlon(ds_sw_cor,mask_holder[slice_buffer], "MASK", "-", "Actual Glacier Mask" ) + + + ds_sw_cor['MASK'] = (('lat','lon'),np.where(ds_sw_cor['MASK'] == 1, ds_sw_cor['MASK'], np.nan)) + if elevation_profile == False: + ds_sw_cor = ds_sw_cor[['sw_dir_cor','MASK']] + print("concat took:", time.time()-now) + + #regrid static ds and merge with sw_dir_cor + if regrid == True: + regrid_no_mask = xe.Regridder(static_ds, ds_coarse[["HGT"]], method="conservative_normed") + regrid = regrid_no_mask(static_ds, ds_coarse[["HGT"]]) + combined = xr.merge([ds_sw_cor, regrid]) + else: + if elevation_profile == True: + combined = ds_sw_cor.copy() + if elevation_profile == False: + combined = xr.merge([ds_sw_cor, static_ds]) + + return combined + + +def run_horayzon_scheme(static_file, file_sw_dir_cor, coarse_static_file=None, + regrid=False, elevation_profile=False, + elev_bandsize=10, elev_stat_file=None): + # ----------------------------------------------------------------------------- + # Prepare data and initialise Terrain class + # ----------------------------------------------------------------------------- + if elevation_profile == True: + print("Routine check. Regrid Option is set to: ", regrid) + print("Setting regrid to False.") + print("Elevation band size is set to: ", elev_bandsize, "m") + regrid = False + + + + # Load high resolution static data + ds = xr.open_dataset(static_file) + elevation = ds["HGT"].values.copy() #else values get overwritten by later line + elevation_original = ds["HGT"].values.copy() + lon = ds["lon"].values + lat = ds["lat"].values + + slice_in, slice_buffer, mask_glacier, mask_glacier_original = compute_and_slice(lat, lon, ds["MASK"].values) + + print("Inner domain size: " + str(elevation[slice_in].shape)) + + #orthometric height (-> height above mean sea level) + elevation_ortho = np.ascontiguousarray(elevation[slice_in]) + + # Compute ellipsoidal heights + elevation += hray.geoid.undulation(lon, lat, geoid="EGM96") # [m] + + offset_0 = slice_in[0].start + offset_1 = slice_in[1].start + + dem_dim_0, dem_dim_1 = elevation.shape + + vert_grid, vec_norm_enu, vec_tilt_enu, trans_ecef2enu,\ + x_enu, y_enu, z_enu, rot_mat_glob2loc = compute_coords(lat, lon, elevation, slice_in) + + # Compute surface enlargement factor + surf_enl_fac = 1.0 / (vec_norm_enu * vec_tilt_enu).sum(axis=2) + print("Surface enlargement factor (min/max): %.3f" % surf_enl_fac.min() + + ", %.3f" % surf_enl_fac.max()) + + # Initialise terrain + mask = np.ones(vec_tilt_enu.shape[:2], dtype=np.uint8) + mask[~mask_glacier] = 0 # mask non-glacier grid cells + + terrain = hray.shadow.Terrain() + #dim_in_0, dim_in_1 = vec_tilt_enu.shape[0], vec_tilt_enu.shape[1] + terrain.initialise(vert_grid, dem_dim_0, dem_dim_1, + offset_0, offset_1, vec_tilt_enu, vec_norm_enu, + surf_enl_fac, mask=mask, elevation=elevation_ortho, + refrac_cor=False) + # -> neglect atmospheric refraction -> effect is weak due to high + # surface elevation and thus low atmospheric surface pressure + + # Load Skyfield data + load.directory = _cfg.paths['static_folder'] + planets = load("de421.bsp") + sun = planets["sun"] + earth = planets["earth"] + loc_or = earth + wgs84.latlon(trans_ecef2enu.lat_or, trans_ecef2enu.lon_or) + # -> position lies on the surface of the ellipsoid by default + + # ----------------------------------------------------------------------------- + # Compute Slope and Aspect + # ----------------------------------------------------------------------------- + # Compute slope (in local ENU coordinates!) + vec_tilt_enu_loc = \ + np.ascontiguousarray(hray.topo_param.slope_vector_meth( + x_enu, y_enu, z_enu, + rot_mat=rot_mat_glob2loc, output_rot=True)[1:-1, 1:-1]) + + # Compute slope angle and aspect (in local ENU coordinates) + slope = np.arccos(vec_tilt_enu_loc[:, :, 2].clip(max=1.0)) + #beware of aspect orientation -> N = 0 in HORAYZON, adjust here + aspect = np.pi / 2.0 - np.arctan2(vec_tilt_enu_loc[:, :, 1], + vec_tilt_enu_loc[:, :, 0]) + aspect[aspect < 0.0] += np.pi * 2.0 # [0.0, 2.0 * np.pi] + + #Create output file for HRZ + static_ds = xr.Dataset() + static_ds.coords['lat'] = lat[slice_buffer[0]] + static_ds.coords['lon'] = lon[slice_buffer[1]] + add_variable_along_latlon(static_ds, elevation_ortho[slice_buffer], "elevation", "m", "Orthometric Height") + add_variable_along_latlon(static_ds, np.rad2deg(slope)[slice_buffer], "slope", "degree", "Slope") + add_variable_along_latlon(static_ds, np.rad2deg(aspect)[slice_buffer], "aspect", "m", "Aspect measured clockwise from North") + add_variable_along_latlon(static_ds, surf_enl_fac[slice_buffer], "surf_enl_fac", "-", "Surface enlargement factor") + + # ----------------------------------------------------------------------------- + # Compute correction factor for direct downward shortwave radiation + # ----------------------------------------------------------------------------- + + # Create time axis + # time in UTC, set timeframe here + time_dt_beg = dt.datetime(2020, 1, 1, 0, 00, tzinfo=dt.timezone.utc) + time_dt_end = dt.datetime(2021, 1, 1, 0, 00, tzinfo=dt.timezone.utc) + dt_step = dt.timedelta(hours=1) + num_ts = int((time_dt_end - time_dt_beg) / dt_step) + ta = [time_dt_beg + dt_step * i for i in range(num_ts)] + + # Add sw dir correction and regrid + comp_time_shadow = [] + sw_dir_cor = np.zeros(vec_tilt_enu.shape[:2], dtype=np.float32) + + ##Load coarse grid + if coarse_static_file != None: + ds_coarse = xr.open_dataset(coarse_static_file) + ds_coarse['mask'] = ds_coarse['MASK'] #prepare for masked regridding + else: + #assign empty var + ds_coarse = None + + ### Build regridder ### + #Create sample dataset to use regridding for + #Create data for first timestep. + ts = load.timescale() + t = ts.from_datetime(ta[0]) + astrometric = loc_or.at(t).observe(sun) + alt, az, d = astrometric.apparent().altaz() + x_sun = d.m * np.cos(alt.radians) * np.sin(az.radians) + y_sun = d.m * np.cos(alt.radians) * np.cos(az.radians) + z_sun = d.m * np.sin(alt.radians) + sun_position = np.array([x_sun, y_sun, z_sun], dtype=np.float32) + + terrain.sw_dir_cor(sun_position, sw_dir_cor) + + ## Construct regridder outside of loop - create empty place holder netcdf + result = xr.Dataset() + result.coords['time'] = [pd.to_datetime(ta[0])] + #ix_latmin-11:ix_latmax+11,ix_lonmin-11:ix_lonmax+11 + result.coords['lat'] = lat[slice_buffer[0]] + result.coords['lon'] = lon[slice_buffer[1]] + sw_holder = np.zeros(shape=[1, lat[slice_buffer[0]].shape[0], lon[slice_buffer[1]].shape[0]]) + sw_holder[0,:,:] = sw_dir_cor[slice_buffer] + mask_crop = mask[slice_buffer] + add_variable_along_timelatlon(result, sw_holder, "sw_dir_cor", "-", "correction factor for direct downward shortwave radiation") + add_variable_along_latlon(result, mask_crop, "mask", "-", "Boolean Glacier Mask") + + #build regridder + if regrid == True: + regrid_mask = xe.Regridder(result, ds_coarse, method="conservative_normed") + + result.close() + del result + + #Iterate over timesteps + + ## Note: this script is slower than the original, because we are creating a dataset + ## for each individual timeframe. This could be optimised when there is no regridding required. + ## To allow for a smooth run of both, we ignore this time constraint, as the script is still comparatively fast + + datasets = [] + for i in range(len(ta)): #loop over timesteps + + t_beg = time.time() + + ts = load.timescale() + t = ts.from_datetime(ta[i]) + print(i) + astrometric = loc_or.at(t).observe(sun) + alt, az, d = astrometric.apparent().altaz() + x_sun = d.m * np.cos(alt.radians) * np.sin(az.radians) + y_sun = d.m * np.cos(alt.radians) * np.cos(az.radians) + z_sun = d.m * np.sin(alt.radians) + sun_position = np.array([x_sun, y_sun, z_sun], dtype=np.float32) + + terrain.sw_dir_cor(sun_position, sw_dir_cor) + + comp_time_shadow.append((time.time() - t_beg)) + + ## first create distributed 2d xarray + result = xr.Dataset() + result.coords['time'] = [pd.to_datetime(ta[i])] + #ix_latmin-11:ix_latmax+11,ix_lonmin-11:ix_lonmax+11 + result.coords['lat'] = lat[slice_buffer[0]] + result.coords['lon'] = lon[slice_buffer[1]] + sw_holder = np.zeros(shape=[1, lat[slice_buffer[0]].shape[0], lon[slice_buffer[1]].shape[0]]) + sw_holder[0,:,:] = sw_dir_cor[slice_buffer] + ## this sets the whole small domain to mask == 1 - might have issues in regridding (considers values outside actual mask) but nvm + mask_crop = mask[slice_buffer] + add_variable_along_timelatlon(result, sw_holder, "sw_dir_cor", "-", "correction factor for direct downward shortwave radiation") + # XESMF regridding requires fieldname "mask" to notice it + add_variable_along_latlon(result, mask_crop, "mask", "-", "Boolean Glacier Mask") + + if elevation_profile == True: + elev_holder = elevation_original[slice_in] #could also use elevation_ortho here + mask_holder = mask_glacier_original[slice_in] + add_variable_along_latlon(result,elev_holder[slice_buffer], "HGT", "m asl", "Surface elevation" ) + ## load actual mask + add_variable_along_latlon(result,mask_holder[slice_buffer], "mask_real", "-", "Actual Glacier Mask" ) + + full_elev_range = result["HGT"].values[result["mask_real"] == 1] + bins = np.arange(np.nanmin(full_elev_range), np.nanmax(full_elev_range)+elev_bandsize, elev_bandsize) + labels = bins[:-1] + elev_bandsize/2 + + placeholder = {} + for var in ["SLOPE","ASPECT","lat","lon"]: + placeholder[var] = calculate_1d_elevationband(ds, "HGT", "MASK", var, elev_bandsize) + + for var in ["sw_dir_cor","mask_real"]: + placeholder[var] = calculate_1d_elevationband(result, "HGT", "mask_real", var, elev_bandsize) + + ## construct the dataframe and xarray dataset + #This is the crudest and most simplest try and here I want to avoid having a 26x26 grid filled with NaNs due to computational time + mask_elev = np.ones_like(placeholder['lat'][:-1]) + ## Suggest all points on glacier + df = pd.DataFrame({'lat':placeholder['lat'][:-1], + 'lon': np.mean(placeholder['lon'][:-1]), #just assign the same value for now for simplicity + 'time': pd.to_datetime(ta[i]), + 'HGT': labels, + 'ASPECT': placeholder['ASPECT'][:-1], + 'SLOPE': placeholder['SLOPE'].data, + 'MASK': mask_elev, + 'N_Points': placeholder["mask_real"].data, + 'sw_dir_cor': placeholder["sw_dir_cor"][0,:].data}) + + #drop the timezone argument from pandas datetime object to ensure fluent conversion into xarray + df['time'] = df['time'].dt.tz_localize(None) + ##sort values by index vars, just in case + df.sort_values(by=["time","lat","lon"], inplace=True) + ## if elevation bins are too small, we will not have a unique index .. manual adjust that + # the latitude information is not really useful anymore if we use HORAYZON, so adjust it + try: + df['lat'] = df['lat'] + df['HGT']*1e-9 + df.set_index(['time','lat','lon'], inplace=True) + print(df) + elev_ds = construct_1d_dataset(df) + except: + df['lat'] = df['lat'] + df['HGT']*1e-8 + df.set_index(['time','lat','lon'], inplace=True) + print(df) + elev_ds = construct_1d_dataset(df) + + if regrid == True: + datasets.append(regrid_mask(result)) + #print("regridding took:", time.time()-now) + else: + if elevation_profile == True: + datasets.append(elev_ds) + elev_ds.close() + del elev_ds + del df + else: + datasets.append(result) + #Close and delete files to free memory + result.close() + del result + + time_tot = np.array(comp_time_shadow).sum() + print("Elapsed time (total / per time step): " + "%.2f" % time_tot + + " , %.2f" % (time_tot / len(ta)) + " s") + + combined = merge_timestep_files(datasets=datasets, regrid=regrid, ds_coarse=ds_coarse, + static_ds=static_ds,elevation_profile=elevation_profile, + mask_glacier_original=mask_glacier_original, + slice_in=slice_in, + slice_buffer=slice_buffer + ) + + #BBox script to crop to minimal extent! + if elevation_profile == True: + combined.to_netcdf(file_sw_dir_cor) + combined[['HGT','ASPECT','SLOPE','MASK','N_Points']].to_netcdf(elev_stat_file) + else: + cropped_combined = combined.where(combined.MASK == 1, drop=True) + cropped_combined.to_netcdf(file_sw_dir_cor) + +#### !! BEWARE: elevation is not the same when using regridding. Files are the same when using 1D approach. + + +def get_user_arguments(parser: argparse.ArgumentParser) -> argparse.Namespace: + """Get user arguments for converting AWS data. + + Args: + parser: An initialised argument parser. + + Returns: + User arguments for conversion. + """ + parser.description = "Create netCDF input file from a .csv file." + parser.prog = __package__ + + # Required arguments + parser.add_argument( + "-s", + "--static", + dest="static_file", + type=str, + metavar="", + required=True, + help="Path to .nc file with static data", + ) + parser.add_argument( + "-o", + "--output", + dest="file_sw_dir_cor", + type=str, + metavar="", + required=True, + help="Path to the resulting netCDF file", + ) + parser.add_argument( + "-c", + "--coarse-static", + dest="coarse_static_file", + type=str, + #const=None, + default=None, + metavar="", + required=False, + help="Path to coarse static file", + ) + parser.add_argument( + "-r", + "--regridding", + type=str2bool, + #const=False, + default=False, + dest="regrid", + required=False, + help="Boolean whether to regrid or not.", + ) + parser.add_argument( + "-e", + "--elevation_prof", + type=str2bool, + #const=False, + default=False, + dest="elevation_profile", + required=False, + help="Boolean whether to calculate 1D elevation bands", + ) + parser.add_argument( + "-es", + "--elevation_size", + type=int, + dest="elev_bandsize", + #const=None, + default=30, + required=False, + help="Integer on the size of the elevation bands in meters", + ) + parser.add_argument( + "-d", + "--elev_data", + dest="elev_stat_file", + type=str, + #const=None, + default=None, + metavar="", + required=False, + help="Left longitude value of the subset", + ) + arguments = parser.parse_args() + + return arguments + + +def load_config(module_name: str) -> tuple: + """Load configuration for module. + + Args: + module_name: Name of this module. + + Returns: + User arguments and configuration parameters. + """ + params = UtilitiesConfig() + arguments = get_user_arguments(params.parser) + params.load(arguments.utilities_path) + params = params.get_config_expansion(name=module_name) + + return arguments, params + + +def main(): + global _args # Yes, it's bad practice + global _cfg + _args, _cfg = load_config(module_name="create_static") + _args, _cfg = load_config(module_name="create_static") + + run_horayzon_scheme( + _args.static_file, + _args.file_sw_dir_cor, + _args.coarse_static_file, + _args.regrid, + _args.elevation_profile, + _args.elev_bandsize, + _args.elev_stat_file + ) + + +if __name__ == "__main__": + main() diff --git a/cosipy/utilities/aws2cosipy/create_HORAYZON_LUT.py b/cosipy/utilities/aws2cosipy/create_HORAYZON_LUT.py deleted file mode 100644 index 630e0eed..00000000 --- a/cosipy/utilities/aws2cosipy/create_HORAYZON_LUT.py +++ /dev/null @@ -1,516 +0,0 @@ -# Description: Compute gridded correction factor for downward direct shortwave -# radiation from given DEM data (~30 m) and -# mask all non-glacier grid cells according to the glacier outline. -# Consider Earth's surface curvature. -# -# Important note: An Earthdata account is required and 'wget' has to be set -# (https://disc.gsfc.nasa.gov/data-access) to download NASADEM -# data successfully. -# -# Source of applied DEM data: https://lpdaac.usgs.gov/products/nasadem_hgtv001/ -# -# Copyright (c) 2022 ETH Zurich, Christian R. Steger -# MIT License - -##WORK IN PROGRESS!## - -### Current TASK: make 1D version runnable - involves probably an averaging of the SW corr. factor per elevation band -### Is the double rgridding a good work? -### OGGM support recent texts to derive the OGGM-based elevation bands etc.? Combine with older Notebook from patrick - -# Load modules -import os -import numpy as np -import subprocess -from netCDF4 import Dataset, date2num -import zipfile -from skyfield.api import load, wgs84 -import time -import fiona -from rasterio.features import rasterize -from rasterio.transform import Affine -from shapely.geometry import shape -import datetime -import datetime as dt -import horayzon as hray -import xarray as xr -import sys -import xesmf as xe -import pandas as pd - -sys.path.append("../..") -from utilities.aws2cosipy.crop_file_to_glacier import crop_file_to_glacier -from utilities.aws2cosipy.aws2cosipyConfig import WRF - -# ----------------------------------------------------------------------------- -# Settings -# ----------------------------------------------------------------------------- - -# set paths -regrid = False #regrid to coarser resolution -elevation_profile = True ## 1D COSIPY -elev_bandsize = 30 #in m -if elevation_profile == True: - print("Routine check. Regrid Option is set to: ", regrid) - print("Setting regrid to False.") - print("Elevation band size is set to: ", elev_bandsize, "m") - regrid = False -ellps = "WGS84" # Earth's surface approximation (sphere, GRS80 or WGS84) -path_out = "../../data/static/HEF/" -file_sw_dir_cor = "LUT_HORAYZON_sw_dir_cor_1D10m.nc" - -static_file = "../../data/static/HEF/HEF_static_raw.nc" #path to high resolution dataset -coarse_static_file = "../../data/static/HEF/HEF_static_agg.nc" #Load coarse grid - -# ---------------------------- -# Some helper functions -# ---------------------------- - -def add_variable_along_timelatlon(ds, var, name, units, long_name): - """ This function adds missing variables to the DATA class """ - if WRF: - ds[name] = (('time','south_north','west_east'), var) - else: - ds[name] = (('time','lat','lon'), var) - ds[name].attrs['units'] = units - ds[name].attrs['long_name'] = long_name - return ds - -def add_variable_along_latlon(ds, var, name, units, long_name): - """ This function self.adds missing variables to the self.DATA class """ - if WRF: - ds[name] = (('south_north','west_east'), var) - else: - ds[name] = (('lat','lon'), var) - ds[name].attrs['units'] = units - ds[name].attrs['long_name'] = long_name - ds[name].encoding['_FillValue'] = -9999 - return ds - -### function to assign attributes to variable -def assign_attrs(ds, name, units, long_name): - ds[name].attrs['units'] = units - ds[name].attrs['long_name'] = long_name - ds[name].attrs['_FillValue'] = -9999 - -#function for mean of circular values strongly inspired by http://webspace.ship.edu/pgmarr/Geo441/Lectures/Lec%201 -def aspect_means(x): - mean_sine = np.nanmean(np.sin(np.radians(x))) - mean_cosine = np.nanmean(np.cos(np.radians(x))) - r = np.sqrt(mean_cosine**2 + mean_sine**2) - cos_mean = mean_cosine/r - sin_mean = mean_sine/r - mean_angle = np.arctan2(sin_mean, cos_mean) - return np.degrees(mean_angle) - - -## ASPECT test raises error.. fix that! And then do some tests to see that it actually works, especially considering the masks and elevation bins -## WhY?? - np.nanmean also doesnt work, so issue maybe with "map" ? - -def calculate_1d_elevationband(xds, elevation_var, mask_var, var_of_interest, elev_bandsize, slice_idx=None): - - ## first mask vals - xds = xds.where(xds[mask_var] == 1, drop=True) - - #test groupby bins - full_elev_range = xds[elevation_var].values[xds[mask_var] == 1] - bins = np.arange(np.nanmin(full_elev_range), np.nanmax(full_elev_range)+elev_bandsize, elev_bandsize) - labels = bins[:-1] + elev_bandsize/2 - - if var_of_interest in ["lat","lon"]: - values = [] - for i in (bins): - sub = xds.where((xds[mask_var] == 1) & (xds[elevation_var] >= i) & (xds[elevation_var] < i + elev_bandsize), drop=True) - values.append(np.nanmean(sub[var_of_interest].values)) - elif var_of_interest == "ASPECT": - elvs = xds[elevation_var].values.flatten()[xds[mask_var].values.flatten() == 1] - aspects = xds[var_of_interest].values.flatten()[xds[mask_var].values.flatten() == 1] - values = [] - for i in (bins): - values.append(aspect_means(aspects[np.logical_and(elvs >= i, elvs < i+elev_bandsize)])) - elif var_of_interest == mask_var: - if slice_idx is None: - values = xds[var_of_interest].groupby_bins(xds[elevation_var], bins, labels=labels, include_lowest=True).sum(skipna=True, min_count=1) - else: - values = xds[var_of_interest][slice_idx].groupby_bins(xds[elevation_var][slice_idx], bins, labels=labels, include_lowest=True).sum(skipna=True, min_count=1) - ## below calculation doesnt work - #elif var_of_interest in ["aspect","ASPECT"]: - # if slice_idx is None: - # values = xds[var_of_interest].groupby_bins(xds[elevation_var], bins, labels=labels, include_lowest=True).map(aspect_means) - # else: - # values = xds[var_of_interest][slice_idx].groupby_bins(xds[elevation_var][slice_idx], bins, labels=labels, include_lowest=True).map(aspect_means) - else: - if slice_idx is None: - values = xds[var_of_interest].groupby_bins(xds[elevation_var], bins, labels=labels, include_lowest=True).mean(skipna=True) - - else: - values = xds[var_of_interest][slice_idx].groupby_bins(xds[elevation_var][slice_idx], bins, labels=labels, include_lowest=True).mean(skipna=True) - - return values - -def construct_1d_dataset(df): - elev_ds = df.to_xarray() - elev_ds.lon.attrs['standard_name'] = 'lon' - elev_ds.lon.attrs['long_name'] = 'longitude' - elev_ds.lon.attrs['units'] = 'Average Lon of elevation bands' - - elev_ds.lat.attrs['standard_name'] = 'lat' - elev_ds.lat.attrs['long_name'] = 'latitude' - elev_ds.lat.attrs['units'] = 'Average Lat of elevation bands' - assign_attrs(elev_ds, 'HGT','meters','Mean of elevation range per bin as meter above sea level') - assign_attrs(elev_ds, 'ASPECT','degrees','Mean Aspect of slope') - assign_attrs(elev_ds, 'SLOPE','degrees','Mean Terrain slope') - assign_attrs(elev_ds, 'MASK','boolean','Glacier mask') - assign_attrs(elev_ds, 'N_Points','count','Number of Points in each bin') - assign_attrs(elev_ds, 'sw_dir_cor','-','Average shortwave radiation correction factor per elevation band') - - return elev_ds - -# ----------------------------------------------------------------------------- -# Prepare data and initialise Terrain class -# ----------------------------------------------------------------------------- - -# Check if output directory exists -if not os.path.isdir(path_out): - os.makedirs(path_out, exist_ok=True) - -# Load high resolution static data -ds = xr.open_dataset(static_file) -elevation = ds["HGT"].values.copy() #else values get overwritten by later line -elevation_original = ds["HGT"].values.copy() -lon = ds["lon"].values -lat = ds["lat"].values - -# Compute indices of inner domain -> needs to encompass everything in range for aggregation -slice_in = (slice(1,lat.shape[0]-1, None), slice(1, lon.shape[0]-1)) - -offset_0 = slice_in[0].start -offset_1 = slice_in[1].start -print("Inner domain size: " + str(elevation[slice_in].shape)) - -#orthometric height (-> height above mean sea level) -elevation_ortho = np.ascontiguousarray(elevation[slice_in]) - -# Compute ellipsoidal heights -elevation += hray.geoid.undulation(lon, lat, geoid="EGM96") # [m] - -# Compute glacier mask -mask_glacier_original = ds["MASK"].values -#set NaNs to zero, relict from create static file -mask_glacier_original[np.isnan(mask_glacier_original)] = 0 -mask_glacier = mask_glacier_original.astype(bool) -mask_glacier = mask_glacier[slice_in] #-1 -1 verywhere - -#mask with buffer for aggregation to lower spatial resolutions -#set +- 11 grid cells to "glacier" to allow ensure regridding -ilist = [] -jlist = [] -## Note that this list is not based on the original shape, see slice_in above -for i in np.arange(0,mask_glacier.shape[0]): - for j in np.arange(0,mask_glacier.shape[1]): - if mask_glacier[i,j] == True: - #print("Grid cell is glacier.") - ilist.append(i) - jlist.append(j) -#create buffer around glacier -ix_latmin = np.min(ilist) -ix_latmax = np.max(ilist) -ix_lonmin = np.min(jlist) -ix_lonmax = np.max(jlist) - -#Watch out that the large domain incorporates the buffer -slice_buffer = (slice(ix_latmin-11,ix_latmax+11), slice(ix_lonmin-11, ix_lonmax+11)) -mask_glacier[slice_buffer] = True - -# Compute ECEF coordinates -x_ecef, y_ecef, z_ecef = hray.transform.lonlat2ecef(*np.meshgrid(lon, lat), - elevation, ellps=ellps) -dem_dim_0, dem_dim_1 = elevation.shape - -# Compute ENU coordinates -trans_ecef2enu = hray.transform.TransformerEcef2enu( - lon_or=lon[int(len(lon) / 2)], lat_or=lat[int(len(lat) / 2)], ellps=ellps) -x_enu, y_enu, z_enu = hray.transform.ecef2enu(x_ecef, y_ecef, z_ecef, - trans_ecef2enu) - -# Compute unit vectors (up and north) in ENU coordinates for inner domain -vec_norm_ecef = hray.direction.surf_norm(*np.meshgrid(lon[slice_in[1]], - lat[slice_in[0]])) -vec_north_ecef = hray.direction.north_dir(x_ecef[slice_in], y_ecef[slice_in], - z_ecef[slice_in], vec_norm_ecef, - ellps=ellps) -del x_ecef, y_ecef, z_ecef -vec_norm_enu = hray.transform.ecef2enu_vector(vec_norm_ecef, trans_ecef2enu) -vec_north_enu = hray.transform.ecef2enu_vector(vec_north_ecef, trans_ecef2enu) -del vec_norm_ecef, vec_north_ecef - -# Merge vertex coordinates and pad geometry buffer -# holds all the data -vert_grid = hray.auxiliary.rearrange_pad_buffer(x_enu, y_enu, z_enu) - -# Compute rotation matrix (global ENU -> local ENU) - -rot_mat_glob2loc = hray.transform.rotation_matrix_glob2loc(vec_north_enu, - vec_norm_enu) - -del vec_north_enu - -# Compute slope (in global ENU coordinates!) -slice_in_a1 = (slice(slice_in[0].start - 1, slice_in[0].stop + 1), - slice(slice_in[1].start - 1, slice_in[1].stop + 1)) - -## Slope vs plain method -> for comparison later -vec_tilt_enu = \ - np.ascontiguousarray(hray.topo_param.slope_vector_meth( - x_enu[slice_in_a1], y_enu[slice_in_a1], z_enu[slice_in_a1], - rot_mat=rot_mat_glob2loc, output_rot=False)[1:-1, 1:-1]) - -# Compute surface enlargement factor -surf_enl_fac = 1.0 / (vec_norm_enu * vec_tilt_enu).sum(axis=2) -print("Surface enlargement factor (min/max): %.3f" % surf_enl_fac.min() - + ", %.3f" % surf_enl_fac.max()) - -# Initialise terrain -mask = np.ones(vec_tilt_enu.shape[:2], dtype=np.uint8) -mask[~mask_glacier] = 0 # mask non-glacier grid cells - -terrain = hray.shadow.Terrain() -dim_in_0, dim_in_1 = vec_tilt_enu.shape[0], vec_tilt_enu.shape[1] -terrain.initialise(vert_grid, dem_dim_0, dem_dim_1, - offset_0, offset_1, vec_tilt_enu, vec_norm_enu, - surf_enl_fac, mask=mask, elevation=elevation_ortho, - refrac_cor=False) -# -> neglect atmospheric refraction -> effect is weak due to high -# surface elevation and thus low atmospheric surface pressure - -# Load Skyfield data -load.directory = path_out -planets = load("de421.bsp") -sun = planets["sun"] -earth = planets["earth"] -loc_or = earth + wgs84.latlon(trans_ecef2enu.lat_or, trans_ecef2enu.lon_or) -# -> position lies on the surface of the ellipsoid by default - -# ----------------------------------------------------------------------------- -# Compute Slope and Aspect -# ----------------------------------------------------------------------------- -# Compute slope (in local ENU coordinates!) -vec_tilt_enu_loc = \ - np.ascontiguousarray(hray.topo_param.slope_vector_meth( - x_enu, y_enu, z_enu, - rot_mat=rot_mat_glob2loc, output_rot=True)[1:-1, 1:-1]) - -# Compute slope angle and aspect (in local ENU coordinates) -slope = np.arccos(vec_tilt_enu_loc[:, :, 2].clip(max=1.0)) -#beware of aspect orientation -> N = 0 in HORAYZON, adjust here -aspect = np.pi / 2.0 - np.arctan2(vec_tilt_enu_loc[:, :, 1], - vec_tilt_enu_loc[:, :, 0]) -aspect[aspect < 0.0] += np.pi * 2.0 # [0.0, 2.0 * np.pi] - -#Create output file for HRZ -static_ds = xr.Dataset() -static_ds.coords['lat'] = lat[slice_buffer[0]] -static_ds.coords['lon'] = lon[slice_buffer[1]] -add_variable_along_latlon(static_ds, elevation_ortho[slice_buffer], "elevation", "m", "Orthometric Height") -add_variable_along_latlon(static_ds, np.rad2deg(slope)[slice_buffer], "slope", "degree", "Slope") -add_variable_along_latlon(static_ds, np.rad2deg(aspect)[slice_buffer], "aspect", "m", "Aspect measured clockwise from North") -add_variable_along_latlon(static_ds, surf_enl_fac[slice_buffer], "surf_enl_fac", "-", "Surface enlargement factor") - -# ----------------------------------------------------------------------------- -# Compute correction factor for direct downward shortwave radiation -# ----------------------------------------------------------------------------- - -# Create time axis -# time in UTC, set timeframe here -time_dt_beg = dt.datetime(2020, 1, 1, 0, 00, tzinfo=dt.timezone.utc) -time_dt_end = dt.datetime(2021, 1, 1, 0, 00, tzinfo=dt.timezone.utc) -dt_step = dt.timedelta(hours=1) -num_ts = int((time_dt_end - time_dt_beg) / dt_step) -ta = [time_dt_beg + dt_step * i for i in range(num_ts)] - -# Add sw dir correction and regrid -comp_time_shadow = [] -sw_dir_cor = np.zeros(vec_tilt_enu.shape[:2], dtype=np.float32) - -##Load coarse grid -ds_coarse = xr.open_dataset(coarse_static_file) -ds_coarse['mask'] = ds_coarse['MASK'] #prepare for masked regridding - -### Build regridder ### -#Create sample dataset to use regridding for -#Create data for first timestep. -ts = load.timescale() -t = ts.from_datetime(ta[0]) -astrometric = loc_or.at(t).observe(sun) -alt, az, d = astrometric.apparent().altaz() -x_sun = d.m * np.cos(alt.radians) * np.sin(az.radians) -y_sun = d.m * np.cos(alt.radians) * np.cos(az.radians) -z_sun = d.m * np.sin(alt.radians) -sun_position = np.array([x_sun, y_sun, z_sun], dtype=np.float32) - -terrain.sw_dir_cor(sun_position, sw_dir_cor) - -## Construct regridder outside of loop - create empty place holder netcdf -result = xr.Dataset() -result.coords['time'] = [pd.to_datetime(ta[0])] -#ix_latmin-11:ix_latmax+11,ix_lonmin-11:ix_lonmax+11 -result.coords['lat'] = lat[slice_buffer[0]] -result.coords['lon'] = lon[slice_buffer[1]] -sw_holder = np.zeros(shape=[1, lat[slice_buffer[0]].shape[0], lon[slice_buffer[1]].shape[0]]) -sw_holder[0,:,:] = sw_dir_cor[slice_buffer] -mask_crop = mask[slice_buffer] -add_variable_along_timelatlon(result, sw_holder, "sw_dir_cor", "-", "correction factor for direct downward shortwave radiation") -add_variable_along_latlon(result, mask_crop, "mask", "-", "Boolean Glacier Mask") - -#build regridder -if regrid == True: - regrid_mask = xe.Regridder(result, ds_coarse, method="conservative_normed") - -result.close() -del result - -## Regridder successfully constructed, close and delete place holder - -#Iterate over timesteps - -datasets = [] -for i in range(len(ta)): #loop over timesteps - - t_beg = time.time() - - ts = load.timescale() - t = ts.from_datetime(ta[i]) - print(i) - astrometric = loc_or.at(t).observe(sun) - alt, az, d = astrometric.apparent().altaz() - x_sun = d.m * np.cos(alt.radians) * np.sin(az.radians) - y_sun = d.m * np.cos(alt.radians) * np.cos(az.radians) - z_sun = d.m * np.sin(alt.radians) - sun_position = np.array([x_sun, y_sun, z_sun], dtype=np.float32) - - terrain.sw_dir_cor(sun_position, sw_dir_cor) - - comp_time_shadow.append((time.time() - t_beg)) - - ## first create distributed 2d xarray - result = xr.Dataset() - result.coords['time'] = [pd.to_datetime(ta[i])] - #ix_latmin-11:ix_latmax+11,ix_lonmin-11:ix_lonmax+11 - result.coords['lat'] = lat[slice_buffer[0]] - result.coords['lon'] = lon[slice_buffer[1]] - sw_holder = np.zeros(shape=[1, lat[slice_buffer[0]].shape[0], lon[slice_buffer[1]].shape[0]]) - sw_holder[0,:,:] = sw_dir_cor[slice_buffer] - ## this sets the whole small domain to mask == 1 - might have issues in regridding (considers values outside actual mask) but nvm - mask_crop = mask[slice_buffer] - add_variable_along_timelatlon(result, sw_holder, "sw_dir_cor", "-", "correction factor for direct downward shortwave radiation") - # XESMF regridding requires fieldname "mask" to notice it - add_variable_along_latlon(result, mask_crop, "mask", "-", "Boolean Glacier Mask") - - if elevation_profile == True: - elev_holder = elevation_original[slice_in] #could also use elevation_ortho here - mask_holder = mask_glacier_original[slice_in] - add_variable_along_latlon(result,elev_holder[slice_buffer], "HGT", "m asl", "Surface elevation" ) - ## load actual mask - add_variable_along_latlon(result,mask_holder[slice_buffer], "mask_real", "-", "Actual Glacier Mask" ) - - full_elev_range = result["HGT"].values[result["mask_real"] == 1] - bins = np.arange(np.nanmin(full_elev_range), np.nanmax(full_elev_range)+elev_bandsize, elev_bandsize) - labels = bins[:-1] + elev_bandsize/2 - - placeholder = {} - for var in ["SLOPE","ASPECT","lat","lon"]: - placeholder[var] = calculate_1d_elevationband(ds, "HGT", "MASK", var, elev_bandsize) - - for var in ["sw_dir_cor","mask_real"]: - placeholder[var] = calculate_1d_elevationband(result, "HGT", "mask_real", var, elev_bandsize) - - ## construct the dataframe and xarray dataset - #This is the crudest and most simplest try and here I want to avoid having a 26x26 grid filled with NaNs due to computational time - mask_elev = np.ones_like(placeholder['lat'][:-1]) - ## Suggest all points on glacier - df = pd.DataFrame({'lat':placeholder['lat'][:-1], - 'lon': np.mean(placeholder['lon'][:-1]), #just assign the same value for now for simplicity - 'time': pd.to_datetime(ta[i]), - 'HGT': labels, - 'ASPECT': placeholder['ASPECT'][:-1], - 'SLOPE': placeholder['SLOPE'].data, - 'MASK': mask_elev, - 'N_Points': placeholder["mask_real"].data, - 'sw_dir_cor': placeholder["sw_dir_cor"][0,:].data}) - - #drop the timezone argument from pandas datetime object to ensure fluent conversion into xarray - df['time'] = df['time'].dt.tz_localize(None) - ##sort values by index vars, just in case - df.sort_values(by=["time","lat","lon"], inplace=True) - df.set_index(['time','lat','lon'], inplace=True) - print(df) - elev_ds = construct_1d_dataset(df) - - now = time.time() - if regrid == True: - datasets.append(regrid_mask(result)) - #print("regridding took:", time.time()-now) - else: - if elevation_profile == True: - datasets.append(elev_ds) - elev_ds.close() - del elev_ds - del df - else: - datasets.append(result) - #Close and delete files to free memory - result.close() - del result - - -#Merge single timestep files -now = time.time() -ds_sw_cor = xr.concat(datasets, dim='time') -ds_sw_cor['time'] = pd.to_datetime(ds_sw_cor['time'].values) -if regrid == True: - ds_sw_cor['MASK'] = ds_coarse['MASK'] #replace with original mask, should have same dimensions -else: - #dont have same dimensions - if elevation_profile == True: - ds_sw_cor['HGT'] = ds_sw_cor['HGT'].isel(time=0) - ds_sw_cor['ASPECT'] = ds_sw_cor['ASPECT'].isel(time=0) - ds_sw_cor['SLOPE'] = ds_sw_cor['SLOPE'].isel(time=0) - ds_sw_cor['MASK'] = ds_sw_cor['MASK'].isel(time=0) - ds_sw_cor['N_Points'] = ds_sw_cor['N_Points'].isel(time=0) - else: - mask_holder = mask_glacier_original[slice_in] - add_variable_along_latlon(ds_sw_cor,mask_holder[slice_buffer], "MASK", "-", "Actual Glacier Mask" ) - - -ds_sw_cor['MASK'] = (('lat','lon'),np.where(ds_sw_cor['MASK'] == 1, ds_sw_cor['MASK'], np.nan)) -if elevation_profile == False: - ds_sw_cor = ds_sw_cor[['sw_dir_cor','MASK']] -print("concat took:", time.time()-now) - -time_tot = np.array(comp_time_shadow).sum() -print("Elapsed time (total / per time step): " + "%.2f" % time_tot - + " , %.2f" % (time_tot / len(ta)) + " s") - -#regrid static ds and merge with sw_dir_cor -if regrid == True: - regrid_no_mask = xe.Regridder(static_ds, ds_coarse[["HGT"]], method="conservative_normed") - regrid = regrid_no_mask(static_ds, ds_coarse[["HGT"]]) - combined = xr.merge([ds_sw_cor, regrid]) -else: - if elevation_profile == True: - combined = ds_sw_cor.copy() - if elevation_profile == False: - combined = xr.merge([ds_sw_cor, static_ds]) - -#BBox script to crop to minimal extent! -if elevation_profile == True: - combined.to_netcdf(path_out+file_sw_dir_cor) - combined[['HGT','ASPECT','SLOPE','MASK','N_Points']].to_netcdf(path_out+"HEF_static_10m_elevbands.nc") -else: - #cropped_combined = crop_file_to_glacier(combined) #results in +2 gridsize somehow - cropped_combined = combined.where(combined.MASK ==1, drop=True) - cropped_combined.to_netcdf(path_out+file_sw_dir_cor) - - -## CURRENT VERSION for the Elevation Profile is taking a longer time to run -## First focus on running at all, then we can check for performance improvements diff --git a/cosipy/utilities/aws2cosipy/crop_file_to_glacier.py b/cosipy/utilities/aws2cosipy/crop_file_to_glacier.py index 995bc6db..e89e739e 100644 --- a/cosipy/utilities/aws2cosipy/crop_file_to_glacier.py +++ b/cosipy/utilities/aws2cosipy/crop_file_to_glacier.py @@ -5,31 +5,34 @@ #np.warnings.filterwarnings('ignore') -sys.path.append('../../') -from utilities.aws2cosipy.aws2cosipyConfig import * +from cosipy.utilities.config_utils import UtilitiesConfig import argparse -def crop_file_to_glacier(ds): - - dic_attrs= {'HGT': ('HGT', 'm', 'Elevation'), - 'ASPECT': ('ASPECT', 'degrees', 'Aspect of slope'), - 'SLOPE': ('SLOPE', 'degrees', 'Terrain slope'), - 'MASK': ('MASK', 'boolean', 'Glacier mask'), - 'T2': ('T2', 'K', 'Temperature at 2 m'), - 'RH2': ('RH2', '%', 'Relative humidity at 2 m'), - 'U2': ('U2', 'm s\u207b\xb9', 'Wind velocity at 2 m'), - 'G': ('G', 'W m\u207b\xb2', 'Incoming shortwave radiation'), - 'PRES': ('PRES', 'hPa', 'Atmospheric Pressure'), - 'N_Points': ('N_Points', 'count','Number of Points in each bin'), - 'RRR': ('RRR', 'mm', 'Total precipitation (liquid+solid)'), - 'SNOWFALL': ('SNOWFALL', 'm', 'Snowfall'), - 'LWin': ('LWin', 'W m\u207b\xb2', 'Incoming longwave radiation'), - 'N': ('N', '%', 'Cloud cover fraction'), - 'sw_dir_cor': ('sw_dir_cor', '-', 'correction factor for direct downward shortwave radiation'), - 'slope': ('slope','degrees', 'Horayzon Slope'), - 'aspect': ('aspect','degrees','Horayzon Aspect measured clockwise from the North'), - 'surf_enl_fac': ('surf_enl_fac','-','Surface enlargement factor'), - 'elevation': ('elevation','m','Orthometric Height')} +_args = None +_cfg = None + +def crop_file_to_glacier(ds, WRF, check_vars): + + dic_attrs= { + "HGT": ("HGT", "m", "Elevation"), + "MASK": ("MASK", "boolean", "Glacier mask"), + "SLOPE": ("SLOPE", "degrees", "Terrain slope"), + "ASPECT": ("ASPECT", "degrees", "Aspect of slope"), + "T2": ("T2", "K", "Air temperature at 2 m"), + "RH2": ("RH2", "%", "Relative humidity at 2 m"), + "U2": ("U2", "m s\u207b\xb9", "Wind velocity at 2 m"), + "PRES": ("PRES", "hPa", "Atmospheric pressure"), + "G": ("G", "W m\u207b\xb2", "Incoming shortwave radiation"), + "RRR": ("RRR", "mm", "Total precipitation"), + "SNOWFALL": ("SNOWFALL", "m", "Snowfall"), + "N": ("N", "-", "Cloud fraction"), + "LWin": ("LWin", "W m\u207b\xb2", "Incoming longwave radiation"), + 'N_Points': ('N_Points', 'count','Number of Points in each bin'), + 'sw_dir_cor': ('sw_dir_cor', '-', 'correction factor for direct downward shortwave radiation'), + 'slope': ('slope','degrees', 'Horayzon Slope'), + 'aspect': ('aspect','degrees','Horayzon Aspect measured clockwise from the North'), + 'surf_enl_fac': ('surf_enl_fac','-','Surface enlargement factor'), + 'elevation': ('elevation','m','Orthometric Height')} dso = ds @@ -39,52 +42,33 @@ def crop_file_to_glacier(ds): for var in list(dso.variables): print(var) arr = bbox_2d_array(dso.MASK.values, dso[var].values, var) - if var in ['lat','latitude','lon','longitude','time','Time']: - if var == 'lat' or var == 'latitude': - dso_mod.coords['lat'] = arr - dso_mod.lat.attrs['standard_name'] = 'lat' - dso_mod.lat.attrs['long_name'] = 'latitude' - dso_mod.lat.attrs['units'] = 'degrees_north' - elif var == 'lon' or var == 'longitude': - dso_mod.coords['lon'] = arr - dso_mod.lon.attrs['standard_name'] = 'lon' - dso_mod.lon.attrs['long_name'] = 'longitude' - dso_mod.lon.attrs['units'] = 'degrees_east' + if var in ['lat','latitude','south_north', 'lon','longitude', 'west_east', 'time','Time']: + if var in ['lat','latitude','south_north']: + dso_mod.coords[var] = arr + dso_mod[var].attrs['standard_name'] = var + dso_mod[var].attrs['long_name'] = 'latitude' + dso_mod[var].attrs['units'] = 'degrees_north' + elif var in ['lon', 'longitude', 'west_east']: + dso_mod.coords[var] = arr + dso_mod[var].attrs['standard_name'] = var + dso_mod[var].attrs['long_name'] = 'longitude' + dso_mod[var].attrs['units'] = 'degrees_east' else: dso_mod.coords['time'] = arr elif var in ['HGT','ASPECT','SLOPE','MASK','N_Points','surf_enl_fac','slope','aspect','elevation']: - add_variable_along_latlon(dso_mod, arr, dic_attrs[var][0], dic_attrs[var][1], dic_attrs[var][2]) + add_variable_along_latlon(dso_mod, arr, dic_attrs[var][0], dic_attrs[var][1], dic_attrs[var][2], WRF) else: - add_variable_along_timelatlon(dso_mod, arr, dic_attrs[var][0], dic_attrs[var][1], dic_attrs[var][2]) + add_variable_along_timelatlon(dso_mod, arr, dic_attrs[var][0], dic_attrs[var][1], dic_attrs[var][2], WRF) #---------------------- # Do some checks #---------------------- print("Performing checks.") - check_for_nan(dso_mod) - - if (T2_var in list(dso_mod.variables)): - check(dso_mod.T2,316.16,223.16) - if (RH2_var in list(dso_mod.variables)): - check(dso_mod.RH2,100.0,0.0) - if (U2_var in list(dso_mod.variables)): - check(dso_mod.U2, 50.0, 0.0) - if (G_var in list(dso_mod.variables)): - check(dso_mod.G,1600.0,0.0) - if (PRES_var in list(dso_mod.variables)): - check(dso_mod.PRES,1080.0,200.0) - - if (RRR_var in list(dso_mod.variables)): - check(dso_mod.RRR,25.0,0.0) - - if (SNOWFALL_var in list(dso_mod.variables)): - check(dso_mod.SNOWFALL, 0.05, 0.0) - - if (LWin_var in list(dso_mod.variables)): - check(dso_mod.LWin, 400, 0.0) - - if (N_var in list(dso_mod.variables)): - check(dso_mod.N, 1.0, 0.0) + check_for_nan(dso_mod, WRF) + + # Check data + if check_vars is True: + check_data(dataset=dso_mod) return dso_mod @@ -102,12 +86,12 @@ def bbox_2d_array(mask, arr, varname): if varname in ['time','Time']: i_min = 0 i_max = None - elif varname in ['lat','latitude']: + elif varname in ['lat','latitude','south_north']: ix = np.where(np.any(mask == 1, axis=1))[0] i_min, i_max = ix[[0, -1]] i_min = i_min #lower bounds included i_max = i_max +1 - elif varname in ['lon','longitude']: + elif varname in ['lon','longitude','west_east']: ix = np.where(np.any(mask == 1, axis=0))[0] i_min, i_max = ix[[0, -1]] i_min = i_min #lower bounds included @@ -131,58 +115,194 @@ def bbox_2d_array(mask, arr, varname): return bbox -def add_variable_along_timelatlon(ds, var, name, units, long_name): - """ This function adds missing variables to the DATA class """ +def add_variable_along_timelatlon(ds, var, name, units, long_name, WRF): + """Add spatiotemporal data to a dataset.""" if WRF: - ds[name] = (('time','south_north','west_east'), var) + ds[name] = (("time", "south_north", "west_east"), var) else: - ds[name] = (('time','lat','lon'), var) - ds[name].attrs['units'] = units - ds[name].attrs['long_name'] = long_name + ds[name] = (("time", "lat", "lon"), var) + ds[name].attrs["units"] = units + ds[name].attrs["long_name"] = long_name return ds -def add_variable_along_latlon(ds, var, name, units, long_name): - """ This function self.adds missing variables to the self.DATA class """ - if WRF: - ds[name] = (('south_north','west_east'), var) +def add_variable_along_latlon(ds, var, name, units, long_name, WRF): + """Add spatial data to a dataset.""" + if WRF: + ds[name] = (("south_north", "west_east"), var) else: - ds[name] = (('lat','lon'), var) - ds[name].attrs['units'] = units - ds[name].attrs['long_name'] = long_name - ds[name].encoding['_FillValue'] = -9999 + ds[name] = (("lat", "lon"), var) + ds[name].attrs["units"] = units + ds[name].attrs["long_name"] = long_name + ds[name].encoding["_FillValue"] = -9999 return ds -def check(field, max, min): - '''Check the validity of the input data ''' - if np.nanmax(field) > max or np.nanmin(field) < min: - print('\n\nWARNING! Please check the data, its seems they are out of a reasonable range %s MAX: %.2f MIN: %.2f \n' % (str.capitalize(field.name), np.nanmax(field), np.nanmin(field))) +def check_data(dataset: xr.Dataset): + """Check data is within physically reasonable bounds.""" + T2_var = _cfg.names['T2_var'] + PRES_var = _cfg.names['PRES_var'] + RH2_var = _cfg.names['RH2_var'] + G_var = _cfg.names['G_var'] + RRR_var = _cfg.names['RRR_var'] + U2_var = _cfg.names['U2_var'] + LWin_var = _cfg.names['LWin_var'] + SNOWFALL_var = _cfg.names['SNOWFALL_var'] + N_var = _cfg.names['N_var'] + + var_list = list(dataset.variables) + + if T2_var in var_list: + check(dataset[T2_var], 316.16, 223.16) + if RH2_var in var_list: + check(dataset[RH2_var], 100.0, 0.0) + if U2_var in var_list: + check(dataset[U2_var], 50.0, 0.0) + if G_var in var_list: + check(dataset[G_var], 1600.0, 0.0) + if PRES_var in var_list: + check(dataset[PRES_var], 1080.0, 200.0) + if RRR_var in var_list: + check(dataset[RRR_var], 25.0, 0.0) + if SNOWFALL_var in var_list: + check(dataset[SNOWFALL_var], 0.05, 0.0) + if LWin_var in var_list: + check(dataset[LWin_var], 400, 0.0) + if N_var in var_list: + check(dataset[N_var], 1.0, 0.0) + + +def check(field, max_bound, min_bound): + """Check the validity of the input data.""" + + if np.nanmax(field) > max_bound or np.nanmin(field) < min_bound: + msg = f"{str.capitalize(field.name)} MAX: {np.nanmax(field):.2f} MIN: {np.nanmin(field):.2f}" + print( + f"\n\nWARNING! Please check the data, it seems they are out of a reasonable range {msg}" + ) + +def raise_nan_error(): + """Raise error if NaNs are in the dataset. + + Raises: + ValueError: There are NaNs in the dataset. + """ + raise ValueError("ERROR! There are NaNs in the dataset.") -def check_for_nan(ds): +def check_for_nan(ds, WRF): if WRF is True: - for y,x in product(range(ds.dims['south_north']),range(ds.dims['west_east'])): + for y, x in product( + range(ds.dims["south_north"]), range(ds.dims["west_east"]) + ): mask = ds.MASK.sel(south_north=y, west_east=x) - if mask==1: - if np.isnan(ds.sel(south_north=y, west_east=x).to_array()).any(): - print('ERROR!!!!!!!!!!! There are NaNs in the dataset') - sys.exit() + if mask == 1: + if np.isnan( + ds.sel(south_north=y, west_east=x).to_array() + ).any(): + raise_nan_error() else: - for y,x in product(range(ds.dims['lat']),range(ds.dims['lon'])): + for y, x in product(range(ds.dims["lat"]), range(ds.dims["lon"])): mask = ds.MASK.isel(lat=y, lon=x) - if mask==1: + if mask == 1: if np.isnan(ds.isel(lat=y, lon=x).to_array()).any(): - print('ERROR!!!!!!!!!!! There are NaNs in the dataset') - sys.exit() + raise_nan_error() -if __name__ == "__main__": - - parser = argparse.ArgumentParser(description="Crop input file to glacier bounding box.") - parser.add_argument('-i', '-input_file', dest='input_file', help='Input file to crop to glacier') - parser.add_argument('-o', '-cosipy_file', dest='cosipy_file', help='Name of resulting COSIPY file') +def empty_user_arguments(parser: argparse.ArgumentParser) -> argparse.Namespace: + """Get user arguments for converting AWS data. + + Args: + parser: An initialised argument parser. - args = parser.parse_args() - print('Read input file %s \n' % (args.input_file)) - ds = xr.open_dataset(args.input_file) - dso_mod = crop_file_to_glacier(ds) + Returns: + User arguments for conversion. + """ + parser.description = "Create static file." + parser.prog = __package__ + + arguments = parser.parse_args() + + return arguments + +def get_user_arguments(parser: argparse.ArgumentParser) -> argparse.Namespace: + """Get user arguments for converting AWS data. + + Args: + parser: An initialised argument parser. + + Returns: + User arguments for conversion. + """ + parser.description = "Create netCDF input file from a .csv file." + parser.prog = __package__ + + # Required arguments + parser.add_argument( + "-i", + "--input", + dest="input_file", + type=str, + metavar="", + required=True, + help="Path to .nc file that needs to be cropped", + ) + parser.add_argument( + "-o", + "--output", + dest="cosipy_file", + type=str, + metavar="", + required=True, + help="Path to the resulting COSIPY netCDF file", + ) + + arguments = parser.parse_args() + + return arguments + +def load_config(module_name: str) -> tuple: + """Load configuration for module. + + Args: + module_name: Name of this module. + + Returns: + User arguments and configuration parameters. + """ + params = UtilitiesConfig() + arguments = get_user_arguments(params.parser) + params.load(arguments.utilities_path) + params = params.get_config_expansion(name=module_name) + + return arguments, params + +def load_config_empty(module_name: str) -> tuple: + """Load configuration for module. + + Args: + module_name: Name of this module. + + Returns: + User arguments and configuration parameters. + """ + params = UtilitiesConfig() + arguments = empty_user_arguments(params.parser) + params.load(arguments.utilities_path) + params = params.get_config_expansion(name=module_name) + + return arguments, params + + +def main(): + global _args # Yes, it's bad practice + global _cfg + _args, _cfg = load_config(module_name="aws2cosipy") + + print('Read input file %s \n' % (_args.input_file)) + ds = xr.open_dataset(_args.input_file) + dso_mod = crop_file_to_glacier(ds, WRF=_cfg.coords["WRF"], check_vars=True) ## write out to file ## print("Writing cropped cosipy file.") - dso_mod.to_netcdf(args.cosipy_file) + dso_mod.to_netcdf(_args.cosipy_file) + +if __name__ == "__main__": + main() + + diff --git a/cosipy/utilities/createStatic/updated_static_file.py b/cosipy/utilities/createStatic/updated_static_file.py index 6fba02c6..5b60c89a 100644 --- a/cosipy/utilities/createStatic/updated_static_file.py +++ b/cosipy/utilities/createStatic/updated_static_file.py @@ -1,45 +1,97 @@ """ - This file reads the DEM of the study site and the shapefile and creates the needed static.nc - TODO: Implement 1D elevation band dataset in here as well in case usage with Moelg/Wohlfahrt scheme is desired +Reads the DEM of the study site and the shapefile and creates a +corresponding static file, ``static.nc``. + +Edit the configuration by supplying a valid .toml file. See the sample +``utilities_config.toml`` for more information. + +Usage: + +From source: +``python -m cosipy.utilities.createStatic.create_static_file [-u ]`` + +Entry point: +``cosipy-create-static [-u ]`` + +Optional arguments: + -u, --u Relative path to utilities' configuration file. """ -import sys +import argparse import os -import xarray as xr -import numpy as np from itertools import product + +import numpy as np import richdem as rd +import xarray as xr import fiona from horayzon.domain import curved_grid -static_folder = '../../data/static/HEF/' +from cosipy.utilities.config_utils import UtilitiesConfig +from cosipy.utilities.aws2cosipy.crop_file_to_glacier import crop_file_to_glacier + +_args = None +_cfg = None + +def check_folder_path(path: str) -> str: + """Check the folder path includes a forward slash.""" + if not path.endswith("/"): + path = f"{path}/" + + return path + -sys.path.append("../..") -from utilities.aws2cosipy.crop_file_to_glacier import crop_file_to_glacier +def check_for_nan(ds,var=None): + for y,x in product(range(ds.dims['lat']),range(ds.dims['lon'])): + mask = ds.MASK.isel(lat=y, lon=x) + if mask==1: + if var is None: + if np.isnan(ds.isel(lat=y, lon=x).to_array()).any(): + raise ValueError("ERROR! There are NaNs in the static fields") + else: + if np.isnan(ds[var].isel(lat=y, lon=x)).any(): + raise ValueError("ERROR! There are NaNs in the static fields") -## Define settings ## -distributed_radiation = True -tile = True -## If distributed radation is not desired, run default static file creation according to values below ## -# to aggregate the DEM to a coarser spatial resolution -aggregate = False -aggregate_degree = '0.00277778' #300m -automatic_domain = True #Do km buffer around glacier or set lat lon box by hand in functions below -crop_file = False #crop to minimum extent -# ref exist: If already have high res. static data set to True and skip calculation -ref_exist = False +def insert_var(ds, var, name, units, long_name): + """Insert variables in dataset""" + ds[name] = (('lat','lon'), var) + ds[name].attrs['units'] = units + ds[name].attrs['long_name'] = long_name + ds[name].attrs['_FillValue'] = -9999 -### input digital elevation model (DEM) -dem_path_tif = static_folder + 'DEM/COPDEM30_HEF_Box.tif' #n30 -### input shape of glacier or study area, e.g. from the Randolph glacier inventory +def get_user_arguments(parser: argparse.ArgumentParser) -> argparse.Namespace: + """Get user arguments for converting AWS data. + + Args: + parser: An initialised argument parser. + + Returns: + User arguments for conversion. + """ + parser.description = "Create static file." + parser.prog = __package__ -shape_path = static_folder + 'Shapefiles/HEF_RGI6.shp' + arguments = parser.parse_args() -### path were the static.nc file is saved -output_path = static_folder + 'HEF_static_raw.nc' -output_path_crop = static_folder + 'HEF_static_raw_crop.nc' -output_path_agg = static_folder + 'HEF_static_agg.nc' + return arguments + + +def load_config(module_name: str) -> tuple: + """Load configuration for module. + + Args: + module_name: name of this module. + + Returns: + User arguments and configuration parameters. + """ + params = UtilitiesConfig() + arguments = get_user_arguments(params.parser) + params.load(arguments.utilities_path) + params = params.get_config_expansion(name=module_name) + + return arguments, params #domain creation assumes WGS84 is valid def domain_creation(shp_path, dist_search, ellps="WGS84"): @@ -60,63 +112,92 @@ def domain_creation(shp_path, dist_search, ellps="WGS84"): return (longitude_upper_left, latitude_upper_left, longitude_lower_right, latitude_lower_right) -def create_static(dem_path_tif=dem_path_tif, shape_path=shape_path, output_path=output_path, output_path_crop=output_path_crop, tile=tile, aggregate=aggregate, aggregate_degree=aggregate_degree, automatic_domain=automatic_domain, crop_file=crop_file, dist_search=20.0): + +def main(): + _, _cfg = load_config(module_name="create_static") + + static_folder = _cfg.paths["static_folder"] + tile = _cfg.coords["tile"] + aggregate = _cfg.coords["aggregate"] + + # input digital elevation model (DEM) + dem_path_tif = f"{static_folder}{_cfg.paths['dem_path']}" + # input shape of glacier or study area, e.g. from the Randolph glacier inventory + shape_path = f"{static_folder}{_cfg.paths['shape_path']}" + # path where the static.nc file is saved + output_path = f"{static_folder}{_cfg.paths['output_file']}" + output_path_crop = f"{static_folder}{_cfg.paths['output_file_crop']}" + + automatic_domain = _cfg.coords["automatic_domain"] + dist_search = _cfg.coords["dist_search"] if automatic_domain: - longitude_upper_left, latitude_upper_left, longitude_lower_right, latitude_lower_right = domain_creation(shape_path, dist_search=dist_search, ellps="WGS84") + longitude_upper_left, latitude_upper_left,\ + longitude_lower_right, latitude_lower_right = domain_creation(shape_path, + dist_search=dist_search, + ellps="WGS84") else: ### to shrink the DEM use the following lat/lon corners print("Please be aware that you switched off the automatic domain creation which requires by-hand adjustment of the borders.") - longitude_upper_left = '90.62' - latitude_upper_left = '30.48' - longitude_lower_right = '90.66' - latitude_lower_right = '30.46' - - ### to aggregate the DEM to a coarser spatial resolution - aggregate_degree = aggregate_degree - - ### intermediate files, will be removed afterwards - dem_path_tif_temp = static_folder + 'DEM_temp.tif' - dem_path_tif_temp2 = static_folder + 'DEM_temp2.tif' - dem_path = static_folder + 'dem.nc' - aspect_path = static_folder + 'aspect.nc' - mask_path = static_folder + 'mask.nc' - slope_path = static_folder + 'slope.nc' - - ### If you do not want to shrink the DEM, comment out the following to three lines + + # to shrink the DEM use the following lat/lon corners + longitude_upper_left = str(_cfg.coords["longitude_upper_left"]) + latitude_upper_left = str(_cfg.coords["latitude_upper_left"]) + longitude_lower_right = str(_cfg.coords["longitude_lower_right"]) + latitude_lower_right = str(_cfg.coords["latitude_lower_right"]) + + # to aggregate the DEM to a coarser spatial resolution + aggregate_degree = str(_cfg.coords["aggregate_degree"]) + + # intermediate files, will be removed afterwards + dem_path_tif_temp = f"{static_folder}DEM_temp.tif" + dem_path_tif_temp2 = f"{static_folder}DEM_temp2.tif" + dem_path = f"{static_folder}dem.nc" + aspect_path = f"{static_folder}aspect.nc" + mask_path = f"{static_folder}mask.nc" + slope_path = f"{static_folder}slope.nc" + if tile: - os.system('gdal_translate -projwin ' + longitude_upper_left + ' ' + latitude_upper_left + ' ' + - longitude_lower_right + ' ' + latitude_lower_right + ' ' + dem_path_tif + ' ' + dem_path_tif_temp) + os.system( + f"gdal_translate -projwin {longitude_upper_left} " + + f"{latitude_upper_left} {longitude_lower_right} " + + f"{latitude_lower_right} {dem_path_tif} {dem_path_tif_temp}" + ) dem_path_tif = dem_path_tif_temp - ### If you do not want to aggregate DEM, comment out the following to two lines if aggregate: - print("Aggregating to ", aggregate_degree) - os.system('gdalwarp -tr ' + aggregate_degree + ' ' + aggregate_degree + ' -r average ' + dem_path_tif + ' ' + dem_path_tif_temp2) + os.system( + f"gdalwarp -tr {aggregate_degree} {aggregate_degree} -r average " + + f"{dem_path_tif} {dem_path_tif_temp2}" + ) dem_path_tif = dem_path_tif_temp2 - ### convert DEM from tif to NetCDF - os.system('gdal_translate -of NETCDF ' + dem_path_tif + ' ' + dem_path) + # convert DEM from tif to NetCDF + os.system(f"gdal_translate -of NETCDF {dem_path_tif} {dem_path}") - ### calculate slope as NetCDF from DEM - os.system('gdaldem slope -of NETCDF ' + dem_path + ' ' + slope_path + ' -s 111120') + # calculate slope as NetCDF from DEM + os.system(f"gdaldem slope -of NETCDF {dem_path} {slope_path} -s 111120") - ### calculate aspect from DEM - #no_data value may differ from file to file, depending on your DEM you may need to change this - aspect = np.flipud(rd.TerrainAttribute(rd.LoadGDAL(dem_path_tif, no_data=-9999), attrib = 'aspect')) + # calculate aspect from DEM + #rd might throw an error when the no_data value is not directly specified + try: + aspect = np.flipud(rd.TerrainAttribute(rd.LoadGDAL(dem_path_tif), attrib = 'aspect')) + except: + aspect = np.flipud(rd.TerrainAttribute(rd.LoadGDAL(dem_path_tif, no_data=-9999), attrib = 'aspect')) - ### calculate mask as NetCDF with DEM and shapefile - os.system('gdalwarp -of NETCDF --config GDALWARP_IGNORE_BAD_CUTLINE YES -cutline ' + shape_path + ' ' + dem_path_tif + ' ' + mask_path) + # calculate mask as NetCDF with DEM and shapefile + os.system( + f"gdalwarp -of NETCDF --config GDALWARP_IGNORE_BAD_CUTLINE YES " + + f"-cutline {shape_path} {dem_path_tif} {mask_path}" + ) - ### open intermediate netcdf files + # open intermediate netcdf files dem = xr.open_dataset(dem_path) mask = xr.open_dataset(mask_path) slope = xr.open_dataset(slope_path) - ### set NaNs in mask to -9999 and elevation within the shape to 1 + # set NaNs in mask to -9999 and elevation within the shape to 1 mask=mask.Band1.values - #some datasets have 0s instead of -9999 why do they occur when doing aggregate? - mask[mask == 0] = -9999 mask[np.isnan(mask)]=-9999 mask[mask>0]=1 print(mask) @@ -133,75 +214,41 @@ def create_static(dem_path_tif=dem_path_tif, shape_path=shape_path, output_path= ds.lat.attrs['long_name'] = 'latitude' ds.lat.attrs['units'] = 'degrees_north' - ### function to insert variables to dataset - def insert_var(ds, var, name, units, long_name): - ds[name] = (('lat','lon'), var) - ds[name].attrs['units'] = units - ds[name].attrs['long_name'] = long_name - ds[name].attrs['_FillValue'] = -9999 - - ### insert needed static variables + # insert needed static variables insert_var(ds, dem.Band1.values,'HGT','meters','meter above sea level') insert_var(ds, aspect,'ASPECT','degrees','Aspect of slope') insert_var(ds, slope.Band1.values,'SLOPE','degrees','Terrain slope') insert_var(ds, mask,'MASK','boolean','Glacier mask') - os.system('rm '+ dem_path + ' ' + mask_path + ' ' + slope_path + ' ' + dem_path_tif_temp + ' '+ dem_path_tif_temp2) - - ### save combined static file, delete intermediate files and print number of glacier grid points - def check_for_nan(ds,var=None): - for y,x in product(range(ds.dims['lat']),range(ds.dims['lon'])): - mask = ds.MASK.isel(lat=y, lon=x) - if mask==1: - if var is None: - if np.isnan(ds.isel(lat=y, lon=x).to_array()).any(): - print('ERROR!!!!!!!!!!! There are NaNs in the static fields') - sys.exit() - else: - if np.isnan(ds[var].isel(lat=y, lon=x)).any(): - print('ERROR!!!!!!!!!!! There are NaNs in the static fields') - sys.exit() + os.system( + f"rm {dem_path} {mask_path} {slope_path} " + + f"{dem_path_tif_temp} {dem_path_tif_temp2}" + ) + + """Save combined static file, delete intermediate files, print + number of glacier grid points.""" check_for_nan(ds) - print(output_path) - if crop_file == True: - ds_crop = crop_file_to_glacier(ds) - if aggregate == True: #saved cropped file only + + # crop file to glacier if desired - but do not set to NaN using xr.where + crop_file = _cfg.coords['crop_file'] + if crop_file: + ds_crop = crop_file_to_glacier(ds, WRF=False, check_vars=False) + if aggregate: #save cropped file only ds_crop.to_netcdf(output_path) - else: #save both files + else: + #uncropped file is needed for e.g., HORAYZON calculation (surrounding terrain) + #cropped file is needed as static data for simulation ds.to_netcdf(output_path) ds_crop.to_netcdf(output_path_crop) ds.close() ds_crop.close() else: ds.to_netcdf(output_path) - ds.close() + ds.close() + print("Study area consists of ", np.nansum(mask[mask==1]), " glacier points") print("Done") -#For some reason, runs are switched. - -if distributed_radiation: - #We need to first get a domain at high resolution and then if aggregate is True: create a second domain with lower resolution to use in later stage - if ref_exist: - print("Skipping calculation of high resolution static file.") - else: - #This needs to have a buffer - create_static(dem_path_tif=dem_path_tif, shape_path=shape_path, output_path=output_path, output_path_crop=output_path_crop, - tile=tile, aggregate=False, aggregate_degree=aggregate_degree, - automatic_domain=True, crop_file=True, dist_search=20.0) - print("\n----------------------------------------") - print("Created high resolution domain for LUTs.") - print("----------------------------------------\n") - #This doesn't really need a buffer so crop it to minimum extent possible next - create_static(dem_path_tif=dem_path_tif, shape_path=shape_path, output_path=output_path_agg, output_path_crop=output_path_crop, - tile=tile, aggregate=True, aggregate_degree=aggregate_degree, - automatic_domain=True, crop_file=True, dist_search=1.0) - print("\n----------------------------------------") - print("Stored aggregated domain for resampling.") - print("----------------------------------------\n") -else: - #Generally less grid cells the better, so crop to minimum extent possible - create_static(dem_path_tif=dem_path_tif, shape_path=shape_path, output_path=output_path, output_path_crop=output_path_crop, - tile=tile, aggregate=False, aggregate_degree=aggregate_degree, - automatic_domain=True, crop_file=True, dist_search=20.0) -#cropfile=False +if __name__ == "__main__": + main() +##old version has NaNs, new version only 0 and 1s \ No newline at end of file