diff --git a/.gitignore b/.gitignore index f80601dd..3f8f01cb 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +.ipynb_checkpoints .idea/ venv/ htmlcov/ diff --git a/cosipy/cpkernel/cosipy_core.py b/cosipy/cpkernel/cosipy_core.py index 1d9cf090..5cfdc108 100644 --- a/cosipy/cpkernel/cosipy_core.py +++ b/cosipy/cpkernel/cosipy_core.py @@ -19,9 +19,10 @@ from cosipy.cpkernel.init import init_snowpack, load_snowpack from cosipy.cpkernel.io import IOClass +from cosipy.utils.options import read_opt -def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_data=None): +def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_data=None, opt_dict=None): """ Cosipy core function, which perform the calculations on one core. Params @@ -49,6 +50,9 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat max_layers = int(DATA.max_layers.values) z = float(DATA.ZLVL.values) + # Replace imported variables with content of the opt_dict. If it's empty + # nothing happens. + read_opt(opt_dict, globals()) # Local variables nt = len(DATA.time.values) #accessing DATA is expensive _RRR = np.full(nt, np.nan) @@ -96,14 +100,14 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat # Initialize snowpack or load restart grid #-------------------------------------------- if GRID_RESTART is None: - GRID = init_snowpack(DATA) + GRID = init_snowpack(DATA, opt_dict) else: GRID = load_snowpack(GRID_RESTART) # Create the local output datasets if not coupled RESTART = None if not WRF_X_CSPY: - IO = IOClass(DATA) + IO = IOClass(DATA, opt_dict) RESTART = IO.create_local_restart_dataset() # hours since the last snowfall (albedo module) @@ -222,12 +226,12 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat #-------------------------------------------- # Calculate albedo and roughness length changes if first layer is snow #-------------------------------------------- - alpha = updateAlbedo(GRID) + alpha = updateAlbedo(GRID, opt_dict) #-------------------------------------------- # Update roughness length #-------------------------------------------- - z0 = updateRoughness(GRID) + z0 = updateRoughness(GRID, opt_dict) #-------------------------------------------- # Surface Energy Balance @@ -237,7 +241,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat # Penetrating SW radiation and subsurface melt if SWnet > 0.0: - subsurface_melt, G_penetrating = penetrating_radiation(GRID, SWnet, dt) + subsurface_melt, G_penetrating = penetrating_radiation(GRID, SWnet, dt, opt_dict) else: subsurface_melt = 0.0 G_penetrating = 0.0 @@ -249,12 +253,14 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat # Find new surface temperature (LW is used from the input file) fun, surface_temperature, lw_radiation_in, lw_radiation_out, sensible_heat_flux, latent_heat_flux, \ ground_heat_flux, rain_heat_flux, rho, Lv, MOL, Cs_t, Cs_q, q0, q2 \ - = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t], sw_radiation_net, U2[t], RAIN, SLOPE, LWin=LWin[t]) + = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t], + sw_radiation_net, U2[t], RAIN, SLOPE, LWin=LWin[t], opt_dict=opt_dict) else: # Find new surface temperature (LW is parametrized using cloud fraction) fun, surface_temperature, lw_radiation_in, lw_radiation_out, sensible_heat_flux, latent_heat_flux, \ ground_heat_flux, rain_heat_flux, rho, Lv, MOL, Cs_t, Cs_q, q0, q2 \ - = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t], sw_radiation_net, U2[t], RAIN, SLOPE, N=N[t]) + = update_surface_temperature(GRID, dt, z, z0, T2[t], RH2[t], PRES[t], + sw_radiation_net, U2[t], RAIN, SLOPE, N=N[t], opt_dict=opt_dict) #-------------------------------------------- # Surface mass fluxes [m w.e.q.] @@ -286,22 +292,22 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat #-------------------------------------------- # Percolation #-------------------------------------------- - Q = percolation(GRID, melt + condensation + RAIN/1000.0 + lwc_from_melted_layers, dt) + Q = percolation(GRID, melt + condensation + RAIN/1000.0 + lwc_from_melted_layers, dt, opt_dict) #-------------------------------------------- # Refreezing #-------------------------------------------- - water_refreezed = refreezing(GRID) + water_refreezed = refreezing(GRID, opt_dict) #-------------------------------------------- # Solve the heat equation #-------------------------------------------- - solveHeatEquation(GRID, dt) + solveHeatEquation(GRID, dt, opt_dict) #-------------------------------------------- # Calculate new density to densification #-------------------------------------------- - densification(GRID, SLOPE, dt) + densification(GRID, SLOPE, dt, opt_dict) #-------------------------------------------- # Calculate mass balance diff --git a/cosipy/cpkernel/init.py b/cosipy/cpkernel/init.py index e6ef4926..779e5666 100644 --- a/cosipy/cpkernel/init.py +++ b/cosipy/cpkernel/init.py @@ -3,10 +3,14 @@ from constants import * from config import * from cosipy.cpkernel.grid import * +from cosipy.utils.options import read_opt -def init_snowpack(DATA): +def init_snowpack(DATA, opt_dict): ''' INITIALIZATION ''' + # Read and set options + read_opt(opt_dict, globals()) + ##-------------------------------------- ## Check for WRF data ##-------------------------------------- diff --git a/cosipy/cpkernel/io.py b/cosipy/cpkernel/io.py index 6c2a12b8..76dabda5 100644 --- a/cosipy/cpkernel/io.py +++ b/cosipy/cpkernel/io.py @@ -13,12 +13,15 @@ from constants import * from config import * import configparser +from cosipy.utils.options import read_opt class IOClass: - def __init__(self, DATA=None): + def __init__(self, DATA=None, opt_dict=None): """ Init IO Class""" + # Read and set options + read_opt(opt_dict, globals()) # Read variable list from file config = configparser.ConfigParser() config.read('./cosipy/output') diff --git a/cosipy/modules/albedo.py b/cosipy/modules/albedo.py index 5348745b..dc23500e 100644 --- a/cosipy/modules/albedo.py +++ b/cosipy/modules/albedo.py @@ -1,9 +1,12 @@ import numpy as np from constants import albedo_method, albedo_fresh_snow, albedo_firn, albedo_ice, \ albedo_mod_snow_aging, albedo_mod_snow_depth, snow_ice_threshold +from cosipy.utils.options import read_opt -def updateAlbedo(GRID): +def updateAlbedo(GRID, opt_dict): """ This methods updates the albedo """ + # Read and set options + read_opt(opt_dict, globals()) albedo_allowed = ['Oerlemans98'] if albedo_method == 'Oerlemans98': alphaMod = method_Oerlemans(GRID) diff --git a/cosipy/modules/densification.py b/cosipy/modules/densification.py index fa7e16a1..151eae7e 100644 --- a/cosipy/modules/densification.py +++ b/cosipy/modules/densification.py @@ -2,13 +2,16 @@ from constants import densification_method, snow_ice_threshold, minimum_snow_layer_height, \ zero_temperature, ice_density from numba import njit +from cosipy.utils.options import read_opt -def densification(GRID,SLOPE,dt): +def densification(GRID,SLOPE,dt, opt_dict=None): """ Densification of the snowpack Args: GRID :: GRID-Structure dt :: integration time """ + # Read and set options + read_opt(opt_dict, globals()) densification_allowed = ['Boone', 'Vionnet', 'empirical', 'constant'] if densification_method == 'Boone': diff --git a/cosipy/modules/evaluation.py b/cosipy/modules/evaluation.py index b641d388..7c0113ce 100644 --- a/cosipy/modules/evaluation.py +++ b/cosipy/modules/evaluation.py @@ -1,11 +1,14 @@ import numpy as np from config import eval_method, obs_type +from cosipy.utils.options import read_opt -def evaluate(stake_names, stake_data, df_): +def evaluate(stake_names, stake_data, df_, opt_dict=None): """ This methods evaluates the simulation with the stake measurements stake_name :: """ - + + # Read and set options + read_opt(opt_dict, globals()) if eval_method == 'rmse': stat = rmse(stake_names, stake_data, df_) else: diff --git a/cosipy/modules/heatEquation.py b/cosipy/modules/heatEquation.py index 0a2481ff..cdece316 100644 --- a/cosipy/modules/heatEquation.py +++ b/cosipy/modules/heatEquation.py @@ -1,8 +1,22 @@ import numpy as np from numba import njit +from cosipy.utils.options import read_opt + +def solveHeatEquation(GRID, dt, opt_dict=None): + + # Read and set options + read_opt(opt_dict, globals()) + + heatEquation_method = 'default' + heatEquation_allowed = ['default'] + + if heatEquation_method == 'default': + heatEquation_default(GRID,dt) + else: + raise ValueError("Heat equation = \"{:s}\" is not allowed, must be one of {:s}".format(heatEquation_method, ", ".join(heatEquation_allowed))) @njit -def solveHeatEquation(GRID, dt): +def heatEquation_default(GRID, dt): """ Solves the heat equation on a non-uniform grid dt :: integration time diff --git a/cosipy/modules/penetratingRadiation.py b/cosipy/modules/penetratingRadiation.py index 6d989519..0dfb89c2 100644 --- a/cosipy/modules/penetratingRadiation.py +++ b/cosipy/modules/penetratingRadiation.py @@ -2,8 +2,13 @@ from constants import penetrating_method, snow_ice_threshold, spec_heat_ice, \ zero_temperature, ice_density, water_density, lat_heat_melting from numba import njit +from cosipy.utils.options import read_opt + +def penetrating_radiation(GRID, SWnet, dt, opt_dict=None): + + # Read and set options + read_opt(opt_dict, globals()) -def penetrating_radiation(GRID, SWnet, dt): penetrating_allowed = ['Bintanja95'] if penetrating_method == 'Bintanja95': subsurface_melt, Si = method_Bintanja(GRID, SWnet, dt) diff --git a/cosipy/modules/percolation.py b/cosipy/modules/percolation.py index b4fd1306..84f64600 100644 --- a/cosipy/modules/percolation.py +++ b/cosipy/modules/percolation.py @@ -1,8 +1,26 @@ import numpy as np from numba import njit +from cosipy.utils.options import read_opt + + +def percolation(GRID, water, dt, opt_dict=None): + + # Read and set options + read_opt(opt_dict, globals()) + + percolation_method = 'bucket' + percolation_allowed = ['bucket'] + + if percolation_method == 'bucket': + Q = bucket(GRID,water,dt) + else: + raise ValueError("Percolation method = \"{:s}\" is not allowed, must be one of {:s}".format(percolation_method, ", ".join(percolation_allowed))) + + return Q + @njit -def percolation(GRID, water, dt): +def bucket(GRID, water, dt): """ Percolation and refreezing of melt water through the snow- and firn pack Args: diff --git a/cosipy/modules/refreezing.py b/cosipy/modules/refreezing.py index e3dae6b6..f3cc9898 100644 --- a/cosipy/modules/refreezing.py +++ b/cosipy/modules/refreezing.py @@ -2,9 +2,25 @@ from constants import zero_temperature, spec_heat_ice, ice_density, \ water_density, lat_heat_melting from numba import njit +from cosipy.utils.options import read_opt -@njit -def refreezing(GRID): +def refreezing(GRID, opt_dict=None): + + # Read and set options + read_opt(opt_dict, globals()) + + refreezing_method = 'default' + refreezing_allowed = ['default'] + + if refreezing_method == 'default': + water_refreezed = refreezing_default(GRID) + else: + raise ValueError("Refreezing method = \"{:s}\" is not allowed, must be one of {:s}".format(refreezing_method, ", ".join(refreezing_allowed))) + + return water_refreezed + +#@njit +def refreezing_default(GRID): # water refreezed water_refreezed = 0.0 diff --git a/cosipy/modules/roughness.py b/cosipy/modules/roughness.py index 23595a22..f54fb1b0 100644 --- a/cosipy/modules/roughness.py +++ b/cosipy/modules/roughness.py @@ -1,8 +1,12 @@ from constants import roughness_method, roughness_fresh_snow, \ roughness_firn, roughness_ice, snow_ice_threshold, \ aging_factor_roughness +from cosipy.utils.options import read_opt -def updateRoughness(GRID): +def updateRoughness(GRID, opt_dict=None): + + # Read and set options + read_opt(opt_dict, globals()) roughness_allowed = ['Moelg12'] if roughness_method == 'Moelg12': diff --git a/cosipy/modules/surfaceTemperature.py b/cosipy/modules/surfaceTemperature.py index 8f8b4f0f..575531ad 100644 --- a/cosipy/modules/surfaceTemperature.py +++ b/cosipy/modules/surfaceTemperature.py @@ -5,9 +5,11 @@ from scipy.optimize import minimize, newton from numba import njit from types import SimpleNamespace +from cosipy.utils.options import read_opt -def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLOPE, LWin=None, N=None): +def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLOPE, LWin=None, + N=None, opt_dict=None): """ This methods updates the surface temperature and returns the surface fluxes Given: @@ -47,6 +49,9 @@ def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLO phi :: Stability correction term [-] """ + # Read and set options + read_opt(opt_dict, globals()) + #Interpolate subsurface temperatures to selected subsurface depths for GHF computation B_Ts = interp_subT(GRID) diff --git a/cosipy/utils/edu_utils.py b/cosipy/utils/edu_utils.py new file mode 100644 index 00000000..fe222b2f --- /dev/null +++ b/cosipy/utils/edu_utils.py @@ -0,0 +1,87 @@ +'''Small utility module to clean up the edu notebooks''' +import os +import sys +import numpy as np +from itertools import product +from dask.distributed import Client, as_completed, progress +from cosipy.cpkernel.cosipy_core import cosipy_core +from cosipy.cpkernel.io import IOClass +from config import northing, easting +from cosipy.utils.options import OPTIONS +import pandas as pd + +def create_IO(opt_dict=None, restart=False): + '''Create io, data, results for the run''' + # Create the essentials. + IO = IOClass(opt_dict=opt_dict) + DATA = IO.create_data_file() + RESULT = IO.create_result_file() + + # If we want the restart as well + if restart: + RESTART = IO.create_restart_file() + + return IO, DATA, RESULT, RESTART + + return IO, DATA, RESULT + + +def run_model(DATA, IO, RESULT, opt_dict=None, RESTART=None, + stake_names=None, df_stakes_data=None): + # List to store the jobs "futures" + futures = [] + # Start initializing the cluster + with Client() as client: + print(f'Model is running, check progress here: Progress is currently disabled') + ny = DATA.dims[northing] + nx = DATA.dims[easting] + # We create the jobs that we then submit to the client. + for y, x in product(range(ny), range(nx)): + # We check if the data is within the glacier. + mask = DATA.MASK.isel(lat=y, lon=x) + if mask == 1: + # Check for faulty values. + if np.isnan(DATA.isel(lat=y, lon=x).to_array()).any(): + print('ERROR!!!!!!!!!!! There are NaNs in the dataset') + sys.exit() + # If all is good, append the job to the list. + futures.append(client.submit(cosipy_core, DATA.isel(lat=y, lon=x), y, x, + stake_names=stake_names, + stake_data=df_stakes_data, opt_dict=opt_dict)) + # Then we can do the jobs. + progress(futures) + # Create numpy arrays which aggregates all local results + IO.create_global_result_arrays() + + + for future in as_completed(futures): + # Here we are unpacking the results from the model run, + # getting ready to save it to our RESULTS dataframe. + indY, indX, local_restart, RAIN, SNOWFALL, LWin, LWout, H, LE, B,\ + QRR, MB, surfMB, Q, SNOWHEIGHT, TOTALHEIGHT, TS, ALBEDO, NLAYERS,\ + ME, intMB, EVAPORATION, SUBLIMATION, CONDENSATION, DEPOSITION,\ + REFREEZE, subM, Z0, surfM, MOL, LAYER_HEIGHT, LAYER_RHO, LAYER_T,\ + LAYER_LWC, LAYER_CC, LAYER_POROSITY, LAYER_ICE_FRACTION,\ + LAYER_IRREDUCIBLE_WATER, LAYER_REFREEZE, stake_names,\ + stat, df_eval = future.result() + + + + IO.copy_local_to_global(indY, indX, RAIN, SNOWFALL, LWin, LWout, H, LE, + B, QRR, MB, surfMB, Q, SNOWHEIGHT, TOTALHEIGHT, + TS, ALBEDO, NLAYERS, ME, intMB, EVAPORATION, + SUBLIMATION, CONDENSATION, DEPOSITION, + REFREEZE, subM, Z0, surfM, MOL, LAYER_HEIGHT, + LAYER_RHO, LAYER_T, LAYER_LWC, LAYER_CC, + LAYER_POROSITY, LAYER_ICE_FRACTION, + LAYER_IRREDUCIBLE_WATER, + LAYER_REFREEZE) + + # Write results to file + IO.write_results_to_file() + print('Finished!') + +def print_options(key=None): + '''Get a pretty print of the options we can change.''' + dict_df = pd.DataFrame.from_dict(OPTIONS, orient='index', columns=['value']) + return dict_df.sort_index() diff --git a/cosipy/utils/options.py b/cosipy/utils/options.py new file mode 100644 index 00000000..35c90d32 --- /dev/null +++ b/cosipy/utils/options.py @@ -0,0 +1,33 @@ +from constants import * +from config import * +# Dict to hold valid options. +OPTIONS = { + 'stability_correction': stability_correction, + 'albedo_method': albedo_method, + 'densification_method': densification_method, + 'penetrating_method': penetrating_method, + 'roughness_method': roughness_method, + 'saturation_water_vapour_method': saturation_water_vapour_method, + 'thermal_conductivity_method': thermal_conductivity_method, + 'sfc_temperature_method': sfc_temperature_method, + 'albedo_fresh_snow': albedo_fresh_snow, + 'albedo_firn': albedo_firn, + 'albedo_ice': albedo_ice, + 'roughness_fresh_snow': roughness_fresh_snow, + 'roughness_ice': roughness_ice, + 'roughness_firn': roughness_firn, + 'time_end': time_end, + 'time_start': time_start, + 'input_netcdf': input_netcdf, + 'output_netcdf': output_netcdf + } + +def read_opt(opt_dict, glob): + """ Reads the opt_dict and overwrites the key-value pairs in glob - the calling function's + globals() dictionary.""" + if opt_dict is not None: + for key in opt_dict: + if key in OPTIONS.keys(): + glob[key] = opt_dict[key] + else: + print(f'ATTENTION: {key} is not a valid option. Default will be used!') diff --git a/data/input/Zhadang/Zhadang_ERA5_2009_dst.nc b/data/input/Zhadang/Zhadang_ERA5_2009_dst.nc new file mode 100644 index 00000000..7cefc90f Binary files /dev/null and b/data/input/Zhadang/Zhadang_ERA5_2009_dst.nc differ diff --git a/notebooks/README.md b/notebooks/README.md new file mode 100644 index 00000000..71d2e039 --- /dev/null +++ b/notebooks/README.md @@ -0,0 +1,3 @@ +# COSIPY tutorials notebooks. +Follow this link to launch the notebooks on MyBinder: +[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/Holmgren825/cosipy/edu_revamp?filepath=notebooks%2Fwelcome.ipynb) diff --git a/notebooks/distributed_run.ipynb b/notebooks/distributed_run.ipynb new file mode 100644 index 00000000..ac8f340f --- /dev/null +++ b/notebooks/distributed_run.ipynb @@ -0,0 +1,271 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "mathematical-travel", + "metadata": {}, + "source": [ + "# A first distributed simulation\n", + "\n", + "Until now, we've only been working with single grid point simulations. This is fine for learning and playing around with different configurations of the model since the runs takes less time. A more typical use case of COSIPY is to run a distributed simulation over the entire surface of the glacier, which gives us a much more detailed view of the glacier." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dutch-motion", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to change the cwd for the ipython session, otherwise COSIPY\n", + "# will look for things in the wrong places.\n", + "import os\n", + "import sys\n", + "# This is not really a good method, if cell is re run we end up in the\n", + "# wrong directory.\n", + "os.chdir('./../')\n", + "sys.path.append(os.getcwd())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "applicable-begin", + "metadata": {}, + "outputs": [], + "source": [ + "from cosipy.utils import edu_utils\n", + "# cfg gives us the NAMELIST\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "macro-dinner", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to tell matplotlib to plot inline\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "several-program", + "metadata": {}, + "source": [ + "Setting up a distributed simulation is fairly easy. All we need is 2D input data (3D with time). This data could either be in the form of gcm/reanalysis output or interpolated from station data. For the tutorials we provide a distributed input file `Zhadang_ERA5_2009_dst.nc` for Zhadang. We can take a look at it" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "speaking-syracuse", + "metadata": {}, + "outputs": [], + "source": [ + "# This is where the data is located.\n", + "input_path = './data/input/'\n", + "with xr.open_dataset(input_path+'Zhadang/Zhadang_ERA5_2009_dst.nc') as ds:\n", + " ds = ds.isel(time=slice(0, -1)).load()\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "coupled-whale", + "metadata": {}, + "source": [ + "The input data contains 91 gridpoints. However, only 17 of these are actually within the glacier boundaries and are used for the run. This input data also covers the whole of 2009, so lets also set up the run for a summer month.\n", + "\n", + "
\n", + "
\n", + " Where is the information about the number gridpoints within the glacier? Click me for a hint\n", + " The points within the glacier are selected based on the MASK variable. You can plot it or \"show\" the data by pressing the storage symbol in the output above. It also possible to count it with the method .count(). Try it below!\n", + "
\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "documented-passage", + "metadata": {}, + "outputs": [], + "source": [ + "# Empty cell for the reader\n" + ] + }, + { + "cell_type": "markdown", + "id": "young-homeless", + "metadata": {}, + "source": [ + "## Running a distributed simulation\n", + "\n", + "Now that we've confirmed that the data is of the right dimension, we can initialize the datasets needed to run the simulation. First we have to add the new input file in the `opt_dict`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "reliable-employee", + "metadata": {}, + "outputs": [], + "source": [ + "# The print_options returns a pandas dataframe so we can index it.\n", + "edu_utils.print_options().loc['input_netcdf']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "handed-washer", + "metadata": {}, + "outputs": [], + "source": [ + "# opt_dict\n", + "opt_dict = dict()\n", + "# We change the input_netcdf\n", + "opt_dict['input_netcdf'] = 'Zhadang/Zhadang_ERA5_2009_dst.nc'\n" + ] + }, + { + "cell_type": "markdown", + "id": "instructional-manor", + "metadata": {}, + "source": [ + "And lets change the start and end dates" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bacterial-joyce", + "metadata": {}, + "outputs": [], + "source": [ + "edu_utils.print_options().loc['time_start']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "electric-chair", + "metadata": {}, + "outputs": [], + "source": [ + "opt_dict['time_start'] = '2009-07-01T00:00'\n", + "opt_dict['time_end'] = '2009-07-31T00:00'" + ] + }, + { + "cell_type": "markdown", + "id": "uniform-advantage", + "metadata": {}, + "source": [ + "With this done we can initialize the datasets needed for running the model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "educational-decade", + "metadata": {}, + "outputs": [], + "source": [ + "IO, DATA, RESULTS = edu_utils.create_IO(opt_dict)" + ] + }, + { + "cell_type": "markdown", + "id": "parallel-anthropology", + "metadata": {}, + "source": [ + "The output above confirms that we have 17 glacier gridpoints in our input data.\n", + "\n", + "## Running the model\n", + "We are now ready to run the model. This is just as simple as in the one dimensional case. In the background the results of each gridpoint is calculated individually, we're basically stacking multiple point simulations next to each other. The `edu_utils.run_model` distributes the work so that gridpoints are executed in parallel. Note however that this still might take some time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "interim-indonesian", + "metadata": {}, + "outputs": [], + "source": [ + "edu_utils.run_model(DATA, IO, RESULTS, opt_dict)" + ] + }, + { + "cell_type": "markdown", + "id": "rural-shadow", + "metadata": {}, + "source": [ + "As before we can take a quick look at the results. First however, we have to reduce the dimensions of the data. A simple way is to select a single time step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ordinary-feeling", + "metadata": {}, + "outputs": [], + "source": [ + "RESULTS.isel(time=1).MB.plot();" + ] + }, + { + "cell_type": "markdown", + "id": "compound-consortium", + "metadata": {}, + "source": [ + "Alternatively we can also reduce one of the dimensions by taking the mean of it, creating a so called Hovmöller diagram" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "internal-interim", + "metadata": {}, + "outputs": [], + "source": [ + "RESULTS.MB.mean(dim='lon').plot();" + ] + }, + { + "cell_type": "markdown", + "id": "facial-madagascar", + "metadata": {}, + "source": [ + "## Next steps\n", + "[Back to overview](welcome.ipynb)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cosipy", + "language": "python", + "name": "cosipy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/first_steps.ipynb b/notebooks/first_steps.ipynb new file mode 100644 index 00000000..09c9c170 --- /dev/null +++ b/notebooks/first_steps.ipynb @@ -0,0 +1,319 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# First steps with COSIPY\n", + "\n", + "This is the first tutorial to get you started with using COSIPY. COSIPY can be run in two ways: Either by running `COSIPY.py` from the command line\n", + "\n", + "`\n", + "$ python COSIPY.py\n", + "`\n", + "\n", + "This will read the model configuration from the files `constants.py` and `config.py`, run the model and save the results to a netcdf file on disk. In this case any customization to the run is done by changing the relevant variable(s) in the previously mentioned files. It is also possible to run COSIPY in an interactive session, or from your own scripts. This is what we are going to go through in these tutorials. In this notebook you'll learn how to:\n", + "\n", + "* Setup a simulation for a single gridpoint\n", + "* Run the simulation and save the results\n", + "* Do simple visualisations of the results using matplotlib\n", + "\n", + "First we have to import the python libraries we are going to use" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Have to change the cwd for the ipython session, otherwise COSIPY\n", + "# will look for things in the wrong places.\n", + "import os\n", + "import sys\n", + "# This is not really a good method, if cell is re run we end up in the\n", + "# wrong directory.\n", + "os.chdir('./../')\n", + "sys.path.append(os.getcwd())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# edu_utils has some quality of life functions\n", + "from cosipy.utils import edu_utils\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Have to tell matplotlib to plot inline\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Create input/output files\n", + "\n", + "First we have to create the input and output datasets for the simulation. This is handled by the `IOClass` of COSIPY, which is not only responsible for creating the [xarray](http://xarray.pydata.org/en/stable/index.html) datasets holding the input and output data, but also for updating these datasets with new data after a run and saving the output to a netcdf file.\n", + "\n", + "Before we initialize the `IO`, it is good to take a look at the variables and constants used by the `IO` module, and also the rest of the model. We can plot the defaults of the most interesting ones using the `print_options` from the `edu_utils`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "edu_utils.print_options()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These are loaded from the files `config.py` and `constants.py` located in the root directory of cosipy. We can change any of these variables to customize our simulation. This can by done changing the value of a variable in the previously mentioned files, however for this to take affect we have to reload the `edu_utils` module. \n", + "\n", + "A much simpler way is to store the variable we want to change in in the `opt_dict`. This is a dictionary which is passed along to the different modules of the model and overwrites the default variable. We can specify new values for any number of variables available in the list above. For instance we can change the variable `time_end` to adjust the temporal extent of our simulation:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Inititalize an empty dictionary.\n", + "opt_dict = dict()\n", + "# Create a key-value pair for time end.\n", + "opt_dict['time_end'] = '2009-01-08T00:00'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are now ready to initialize the `IOClass` and the datasets holding the input and output data for our simulation. For this we have a utility function `create_IO`, which takes the `opt_dict` and returns an instance of the `IOClass` and the `DATA` and `RESULTS` datasets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "IO, DATA, RESULTS = edu_utils.create_IO(opt_dict)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The output above gives us some helpful insight about the data. We can see that the data covers the period 2009-01-01 to 2009-01-31, and that we have specified the integration time to 2009-01-01 to 2009-01-08 by adding a new value to `time_end` in the `opt_dict`. We also see that we have one gridpoint which we will perform the calculations on.\n", + "\n", + "It also shows some quality control on the input data, more on this later.\n", + "\n", + "
\n", + " Note that we don't have to use the opt_dict as we did above - if we don't, the simulation will run using the default parameters from config.py and constants.py.\n", + "\n", + "Also note that the IOClass along with the DATA and RESULTS datasets has to be re-initliazed after a variable has been changed/added to the opt_dict.\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Running COSIPY\n", + "Now we are ready to do a simple point simulation with COSIPY. This example is run on the Zhadang glacier. We can change this by changing the `input_netcdf` via the `opt_dict` to point at the files we want to use. In these tutorials data for Zhadang and Hintereisferner are included.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then run the model by simply calling the `run_model` function from the `edu_utils`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + " Note that this can take some time.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "edu_utils.run_model(DATA, IO, RESULTS, opt_dict)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Thats it! Now we are ready to inspect the results." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualizing the results\n", + "To get a quick overview of the variables in our dataframe, we can simply call it " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "RESULTS" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that our data is stored in an xarray dataset - a very handy format for interacting with netcdf files in python. To check which data variables are available, one can simply click on the label \"Data Variables\", which will expand a list of the variables.\n", + "\n", + "We can quckly plot a variable, the surface mass balance (`surfMB`) for instance:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "RESULTS.surfMB.plot();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the use of xarray-datasets for storing the data makes it possible to levarage the built in plotting functionality of xarray and quickly visualize the output.\n", + "\n", + "For a slightly more advanced plot, one can still make use of matplotlib" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, (ax1, ax2) = plt.subplots(figsize=(7,8), nrows=2, sharex=True,\n", + " gridspec_kw = {'height_ratios':[2,1]})\n", + "RESULTS.surfMB.plot(ax=ax1)\n", + "RESULTS.TS.plot(ax=ax2, color='C1')\n", + "# Easier to make use of xarray plotting capabilites and then remove unwanted lables.\n", + "ax1.set_xlabel('')\n", + "ax2.set_title('')\n", + "# Get the time span of our dataset to use for title\n", + "time_span = RESULTS.time[0].values, RESULTS.time[-1].values\n", + "time_span = np.datetime_as_string(time_span, unit='m')\n", + "# Set the title of the plot.\n", + "ax1.set_title(f'Surface mass balance and Snowfall at Zhadang \\n between {time_span[0]} and' +\\\n", + " f' {time_span[1]}')\n", + "# Fixes the spacing between parts of the plots.\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Notice the scale of the surface mass balance.** In the plot above, we see the surface mass blance (blue) and the 2-meter air temperature (orange). Even though the air temperautre stays below 0°C (273.15K) for the entire period, the mass balance is still negative during the day. We will investigate this is the following section." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + " Question: Can you think of a variable to plot as a way to explain the positive surface mass balance during the 6th of January? Click me for a hint\n", + " Try plotting the snowfall, i.e. RESULTS.SNOWFALL.plot().\n", + "
\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Write your code here\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A look at what is causing the melt\n", + "\n", + "Since the temperature never reaches above zero, melt is not causing the small loss of mass during the days seen in the plot above. Instead it is sublimation, the process when a substance transitions from solid to gas without passing through the liquid state. We can plot the sublimation quickly to confirm this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "RESULTS.SUBLIMATION.plot();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Next steps\n", + "[Back to overview](welcome.ipynb)\n", + "\n", + "[Sensitivity studies with COSIPY](sensitivity_study.ipynb)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cosipy", + "language": "python", + "name": "cosipy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.13" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/seb_deep_dive.ipynb b/notebooks/seb_deep_dive.ipynb new file mode 100644 index 00000000..cbc29627 --- /dev/null +++ b/notebooks/seb_deep_dive.ipynb @@ -0,0 +1,429 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "curious-agriculture", + "metadata": {}, + "source": [ + "# A deep dive into the surface mass balance computation in COSIPY" + ] + }, + { + "cell_type": "markdown", + "id": "wooden-nickel", + "metadata": {}, + "source": [ + "In the previous notebooks we've seen how changing different parameters of the model affects the calculation of surface mass balance on our glacier. In order to strengthen our understanding somewhat, in this notebook we will go through how COSIPY is calculating the surface mass balance during a simulation. Lets start from the top.\n", + "\n", + "First we run a simple simulation to have some data to look at." + ] + }, + { + "cell_type": "markdown", + "id": "surprised-hughes", + "metadata": {}, + "source": [ + "**The standard imports**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "sticky-berkeley", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to change the cwd for the ipython session, otherwise COSIPY\n", + "# will look for things in the wrong places.\n", + "import os\n", + "import sys\n", + "# This is not really a good method, if cell is re run we end up in the\n", + "# wrong directory.\n", + "os.chdir('./../')\n", + "sys.path.append(os.getcwd())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "hawaiian-cosmetic", + "metadata": {}, + "outputs": [], + "source": [ + "from cosipy.utils import edu_utils\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "demonstrated-dairy", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to tell matplotlib to plot inline\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "transparent-consolidation", + "metadata": {}, + "source": [ + "Then we initiate the IO and datasets and run the simulation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "second-wholesale", + "metadata": {}, + "outputs": [], + "source": [ + "# IO and datasetes\n", + "IO, DATA, RESULTS = edu_utils.create_IO()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "capital-fraud", + "metadata": {}, + "outputs": [], + "source": [ + "edu_utils.run_model(DATA=DATA, IO=IO, RESULT=RESULTS)" + ] + }, + { + "cell_type": "markdown", + "id": "enabling-signature", + "metadata": {}, + "source": [ + "***\n", + "## Decomposing the mass balance\n", + "The surface mass balance, which is measured in meter water equivalent (m w.e.), is the result of adding the snowfall ($SF$) and deposition ($D$) together and removing the melt ($M$), the sublimation ($s$) and the evaporation/condensation ($e$)\n", + "\n", + "$$\n", + "MB_{surf} = SF + D - M - s - e.\n", + "$$\n", + "\n", + "Note that the evaporation or condensation depends on the sign, it can have either a positive or negative impact on the mass balance We can plot these variables from our `RESULTS` dataset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "crude-swing", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 5))\n", + "mb_vars = ['SNOWFALL', 'surfM', 'SUBLIMATION', 'DEPOSITION', 'EVAPORATION',\n", + " 'CONDENSATION']\n", + "for var in mb_vars:\n", + " RESULTS[var].plot(ax=ax, label=var)\n", + "ax.set_ylabel('m w.e.') \n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "dominican-approach", + "metadata": {}, + "source": [ + "For our simulated period the plot shows us that sublimation and snowfall make up the main part of the changing mass balance. But there is also a small contribution from deposition. Others are negligible.\n", + "\n", + "### Snowfall\n", + "The **snowfall** can be calculated in two ways. If it is a direct product of the weather station/gcm data, it is only multiplied by a constant (`mult_factor_RRR`) and then used directly in the surface mass balance calculation. \n", + "\n", + "In the other case, when only the rain rates are provided, **snowfall** is derived using a logistic transfer function\n", + "\n", + "$$\n", + "\\mathrm{SNOWFALL} = \\frac{RRR}{1000.0} \\frac{\\rho_{ice}}{\\rho_{fresh\\space snow}} \\frac{-tanh((T2 - t_0) - c_1) \\cdot c_2) + 1}{2}\n", + "$$\n", + "\n", + "where $c_1$ and $c_2$ are the center and spread transfer functions. This smoothly scales the solid precipitation from 100% at $t_0$ and 0% at 2 °C. The density of fresh snow is calculated as a function of the air temperature and the wind velocity." + ] + }, + { + "cell_type": "markdown", + "id": "suburban-judge", + "metadata": {}, + "source": [ + "### Melt\n", + "\n", + "The surface melt is calculated from the energy available for melt ($q_m$) and the latent heat required for melting ice ($L_s$)\n", + "\n", + "$$\n", + "\\mathrm{M} = \\frac{q_m \\cdot dt}{L_s \\cdot 1000}.\n", + "$$\n", + "\n", + "We multiply with dt to get the total energy available during the time step, and divide with $1000$ to convert the result to m w.e.\n", + "\n", + "The energy available for melt is the sum of the radiative fluxes:\n", + "- Net short wave radiation\n", + "- Net long wave radiation\n", + "- Ground heat flux\n", + "- Rain heat flux \n", + "- Sensible heat flux\n", + "- Latent heat flux\n", + "\n", + "These are available in the `RESULTS` as well so we can take a look at them" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "indonesian-miracle", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 5))\n", + "(RESULTS['LWin'] + RESULTS['LWout']).plot(ax=ax, label='LW net')\n", + "# Short wave net is actually not in the Results but we can calculate it\n", + "(RESULTS['G'] * (1 - RESULTS['ALBEDO'])).plot(ax=ax, label='SW net')\n", + "RESULTS['H'].plot(ax=ax, label='Sensible heat flux')\n", + "RESULTS['LE'].plot(ax=ax, label='Latent heat flux')\n", + "RESULTS['B'].plot(ax=ax, label='Ground heat flux')\n", + "RESULTS['QRR'].plot(ax=ax, label='Rain heat flux')\n", + "ax.set_ylabel('Energy [W m$^{-2}$]')\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "distinct-butterfly", + "metadata": {}, + "source": [ + "This plot contains a lot of information but there are some key points\n", + "- During the day, the **net short wave** radiation (orange), which act to heat the surface of the glacier is largely compensated by the negative **ground heat flux**. This means that energy is transported into the glacier, slowly heating it.\n", + "- During the night, the sign of the **ground heat flux** is reversed when the glacier release energy in the form of outgoing **long wave radiation**.\n", + "\n", + "From the look of it, the net energy of the variables above should be close to zero. Let's take a look" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "explicit-vision", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 5))\n", + "RESULTS['ME'].plot(ax=ax);" + ] + }, + { + "cell_type": "markdown", + "id": "awful-synthetic", + "metadata": {}, + "source": [ + "With this in mind we can dissect the radiative fluxes even further. \n", + "\n", + "#### Net shortwave radiation and albedo\n", + "The **net shortwave radiation** is defined as\n", + "$$\n", + "q_{sw} = (1-\\alpha) \\cdot q_G\n", + "$$\n", + "\n", + "where $q_G$ is the incoming shortwave radiation and $\\alpha$ is the **albedo**. $q_G$ is given by the input data, while the evolution of the albedo is parametrized according to the approach from Oerlemans and Knap ([1998](https://www.cambridge.org/core/journals/journal-of-glaciology/article/1-year-record-of-global-radiation-and-albedo-in-the-ablation-zone-of-morteratschgletscher-switzerland/F30150EB1AD9D62A0C007C22FFAE7B6A)). In it, the snow albedo is a function of the time since the last snowfall ($s$) and an albedo timescale ($\\tau^*$) that specifies the speed of the degradation of the albedo from that of fresh snow ($\\alpha_s$) to that of firn ($\\alpha_f$)\n", + "\n", + "$$\n", + "\\alpha_{snow} = \\alpha_f + (\\alpha_s - \\alpha_f)\\space exp \\left(\\frac{s}{\\tau^*}\\right).\n", + "$$\n", + "\n", + "As the thickness of the snowpack ($d$) decreases, the surface albedo should approach that of the albedo of ice ($\\alpha_i$). This is done by introducing a characteristic snow depth scale ($d^*$) to employ what is known as e-folding. With this, the full surface albedo can be written as \n", + "\n", + "$$\n", + "\\alpha = \\alpha_s + (\\alpha_i - \\alpha_s)\\space exp \\left(\\frac{-d}{d^*}\\right).\n", + "$$\n", + "\n", + "The surface albedo is reset to that of fresh snow when new snow accumulation exceeds a certain threshold.\n", + "\n", + "The fraction of the net shortwave radiation that penetrates the glacier decays with the depth ($z$) from the surface according to\n", + "\n", + "$$\n", + "Q(z, t) = \\lambda_r q_{sw} \\space exp(-z \\beta)\n", + "$$\n", + "\n", + "where $\\lambda_r$ is the fraction of absorbed radiation and $\\beta$ the extinction coefficient (both depending on if the integrated depth is snow or ice). \n", + "\n", + "
\n", + " Some of the constants in the equations above are possible to configure using the opt_dict. Remember that we changed one of them, the albedo of fresh snow, in the sensitivity study. \n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "waiting-duncan", + "metadata": {}, + "source": [ + "#### Net longwave radiation\n", + "\n", + "The **net longwave radiation** can be handled in two different ways:\n", + "If the input data contains the incoming longwave radiation, the net is calculated as\n", + "\n", + "$$\n", + "q_{lw} = q_{lw_{in}} - \\epsilon_s \\sigma T_0^4,\n", + "$$\n", + "\n", + "where $\\epsilon_s$ is the surface emissivity, in our case constant $\\sim 1$.\n", + "\n", + "If the input data doesn't contain the incoming longwave radiation, this flux will be parametrized using the Stefan-Boltzman law with the means of the air temperature ($T_{zt}$) and the atmospheric emissivity\n", + "\n", + "$$\n", + "\\epsilon_a = \\epsilon_{cs}(1-N^2) + \\epsilon_{cl}N^2,\n", + "$$\n", + "\n", + "where $N$ is the cloud cover fraction and $\\epsilon_{cs}$ and $\\epsilon_{cl}$ are the clear-sky and cloud emissivity. The cloud emissivity is kept constant while the clear-sky emissivity is given by\n", + "\n", + "$$\\epsilon_{cs} = 0.23 + 0.433(e_{zt}/T_{zt})^{1/8}.$$\n", + "\n", + "where $e_{zt}$ is the water vapour pressure.\n", + "\n", + "In our case, the input data does not contain any incoming longwave radiation and we thus rely on the cloud cover fraction and the air temperature for the parameterization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1318ce87-88dc-49ff-b3ec-41795acf35a6", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 5))\n", + "lns1 = (RESULTS['LWin'] + RESULTS['LWout']).plot(ax=ax, label='LW net', c='C0')\n", + "ax2 = ax.twinx()\n", + "lns2 = DATA.N.plot(ax=ax2, c='C1', label='N')\n", + "ax2.set_title('')\n", + "lns = lns1+lns2\n", + "labs = [l.get_label() for l in lns]\n", + "ax.set_ylabel('Energy [W m$^{-2}$]')\n", + "ax.legend(lns, labs);" + ] + }, + { + "cell_type": "markdown", + "id": "duplicate-cooperative", + "metadata": {}, + "source": [ + "#### Sensible and latent heat flux\n", + "\n", + "The sensible and latent heat flux ($q_{sh}$ & $q_{lh}$) are parametrized using a bulk approach of the flux-gradient similarity. This is done since the model only uses meteorological variables from one height. By using the Stanton number ($C_H$) and Dalton number ($C_E$) the turbulent fluxes can be calculated as\n", + "$$\n", + "q_{sh} = \\rho_a c_p C_H u_{z_v} (T_{z_t} - T_0)\n", + "$$\n", + "and\n", + "$$\n", + "q_{lh} = \\rho_a L_v C_E u_{z_v} (q_{z_q} - q_0)\n", + "$$\n", + "with the following constants \n", + "\n", + "|Symbol | Details|\n", + "|:---|:---|\n", + "|$\\rho_a$| Air density|\n", + "|$c_p$| Specific heat of air|\n", + "|$L_v$| Latent heat of vaporisation, (replaced by latent
heat of sublimation $L_s$ if $T_0 < T_m$)|\n", + "|$u_{z_v}$| Wind velocity at height $z_t$|\n", + "|$T_{z_t}$| Temperature at height $z_t$|\n", + "|$q_{z_t}$| Mixing ratio at height $z_t$|\n", + "|$q_{0}$| Mixing ratio at the surface|" + ] + }, + { + "cell_type": "markdown", + "id": "cd63d40e-d6e3-48aa-983e-15b0596f0f96", + "metadata": {}, + "source": [ + "The bulk coefficients ($C_H$ & $C_E$) are independent of the wind speed and only depend on the stability of the atmosphere and the roughness of the surface, $z_{0_v}$. For snow the roughness length is a function of time, which increases linearly from that of fresh snow to that of firn. For ice the roughness length is constant.\n", + "\n", + "\n", + "
\n", + "
\n", + " With your knowledge from the previous notebooks, can you design an experiment to investigate the dependence of the roughness length on the calculations of the turbulent fluxes?\n", + " Click me for a hint!\n", + " Begin by taking a look at the options (edu_utils.print_options()). Then use the appropriate variables for a sensitivity study, much like in the previous notebooks.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dcf4eb57-be23-40a6-99a1-1cbab3fd5da9", + "metadata": {}, + "outputs": [], + "source": [ + "# Write your code here.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26bd36f7-a5d4-4666-a2f2-a7355e4ab0e2", + "metadata": {}, + "outputs": [], + "source": [ + "# Write your code here.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15b4bcf0-968f-4bdd-b240-dd99fd42601d", + "metadata": {}, + "outputs": [], + "source": [ + "# Write your code here. You can add more cells if needed.\n" + ] + }, + { + "cell_type": "markdown", + "id": "overall-dryer", + "metadata": {}, + "source": [ + "#### Ground heat flux" + ] + }, + { + "cell_type": "markdown", + "id": "suffering-variable", + "metadata": {}, + "source": [ + "### Sublimation\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "acting-furniture", + "metadata": {}, + "source": [ + "## Next steps\n", + "[Back to overview](welcome.ipynb)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/sensitivity_study.ipynb b/notebooks/sensitivity_study.ipynb new file mode 100644 index 00000000..a13cf31a --- /dev/null +++ b/notebooks/sensitivity_study.ipynb @@ -0,0 +1,396 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "selective-rapid", + "metadata": {}, + "source": [ + "# Sensitivity studies using COSIPY\n", + "\n", + "In this tutorial we will go through how to do a sensitivity study of the surface mass balance of a glacier using COSIPY. For a sensitivity study we want to check how much changing a constant affects the calculation of a certain variable. We can do this on essentially any variable, but in our case we will focus on the surface mass balance of the glacier.\n", + "\n", + "Just as in [First steps in COSIPY](first_steps.ipynb) we begin with importing some useful libraries and the cosipy functions we need." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "oriented-mining", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to change the cwd for the ipython session, otherwise COSIPY\n", + "# will look for things in the wrong places.\n", + "import os\n", + "import sys\n", + "# This is not really a good method, if cell is re run we end up in the\n", + "# wrong directory.\n", + "os.chdir('./../')\n", + "sys.path.append(os.getcwd())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "heavy-duration", + "metadata": {}, + "outputs": [], + "source": [ + "from cosipy.utils import edu_utils\n", + "# cfg gives us the NAMELIST\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "humanitarian-equality", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to tell matplotlib to plot inline\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "impaired-dragon", + "metadata": {}, + "source": [ + "## Choosing the constant to investigate\n", + "\n", + "Before we can continue with setting up the simulations, we have to choose the constant we want to change. As you might remember, the configurable constants used by COSIPY can be displayed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "mathematical-announcement", + "metadata": {}, + "outputs": [], + "source": [ + "edu_utils.print_options()" + ] + }, + { + "cell_type": "markdown", + "id": "opening-opposition", + "metadata": {}, + "source": [ + "This contains a lot of different things and can be a bit overwhelming at first. It contains variables used in the configuration of where to store files, constants used in some of the sub-modules, and the parameterizations used.\n", + "\n", + "A good guess for a variable affecting the surface mass balance in the winter is the albedo of fresh snow" + ] + }, + { + "cell_type": "markdown", + "id": "considered-monster", + "metadata": {}, + "source": [ + "Let's say we want to change this 10% in either direction. Since the model is overwriting the constants from `opt_dict`, the plan is to change the albedo of fresh snow in there before we use it to initialize the model." + ] + }, + { + "cell_type": "markdown", + "id": "significant-timer", + "metadata": {}, + "source": [ + "### Initializing the opt_dicts and data files\n", + "\n", + "For this experiment we want to do three runs, one default and two with nudged albedo values. We begin by initializing two instances of the `opt_dict`. Remember that we don't need an `opt_dict` for the default run. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "verified-commission", + "metadata": {}, + "outputs": [], + "source": [ + "# First we fetch the default albedo.\n", + "default_albedo = edu_utils.OPTIONS['albedo_fresh_snow']\n", + "# Up 10%\n", + "opt_dict_up = {'albedo_fresh_snow': default_albedo * 1.1}\n", + "# Down 10%\n", + "opt_dict_dn = {'albedo_fresh_snow': default_albedo * 0.9} " + ] + }, + { + "cell_type": "markdown", + "id": "apart-collect", + "metadata": {}, + "source": [ + "Did it work?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "sapphire-olive", + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Up: {opt_dict_up['albedo_fresh_snow']:.3f}\\nDefault: \\\n", + "{default_albedo:.3f}\\nDown: \\\n", + "{opt_dict_dn['albedo_fresh_snow']:.3f}\\n\")" + ] + }, + { + "cell_type": "markdown", + "id": "regulation-reservation", + "metadata": {}, + "source": [ + "
\n", + " Question: Can you figure out what to change in the cells above to alter the albedo by 5% instead of 10%?\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "optimum-championship", + "metadata": {}, + "source": [ + "We can now go on to creating the files used for storing the data and results of the simulation. Just as with the namlist we need one set of instances for each experiment we want to run, in this case three. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "durable-premiere", + "metadata": {}, + "outputs": [], + "source": [ + "# Default\n", + "IO_def, DATA_def, RESULTS_def = edu_utils.create_IO()\n", + "# Up\n", + "IO_up, DATA_up, RESULTS_up = edu_utils.create_IO(opt_dict_up)\n", + "# Down\n", + "IO_dn, DATA_dn, RESULTS_dn = edu_utils.create_IO(opt_dict_dn)" + ] + }, + { + "cell_type": "markdown", + "id": "dependent-tanzania", + "metadata": {}, + "source": [ + "## Running multiple simulations\n", + "\n", + "We have now prepared the namelists and data/results files for the three different runs we want to do. The runs can initilaized as it was done in the [First steps in COSIPY](first_steps.ipynb) notebook, and changing the namelist and data/results variables manually for each run. But to save us writing the same code multiple times, and to give some inspiration how one might go about running a higher number of simulations, we will now do it in a loop. First we have to put all we need into a list" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "lesbian-circumstances", + "metadata": {}, + "outputs": [], + "source": [ + "# List of lists with our experiments\n", + "exp_list = [[DATA_def, IO_def, RESULTS_def, None],\n", + " [DATA_dn, IO_dn, RESULTS_dn, opt_dict_dn],\n", + " [DATA_up, IO_up, RESULTS_up, opt_dict_up]\n", + " ]" + ] + }, + { + "cell_type": "markdown", + "id": "lucky-transfer", + "metadata": {}, + "source": [ + "We can index this list to get the results of experiment 2 as `exp_list[1][2]`." + ] + }, + { + "cell_type": "markdown", + "id": "aerial-hydrogen", + "metadata": {}, + "source": [ + "
\n", + " Remember that python is zero-indexed.\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "identified-geneva", + "metadata": {}, + "source": [ + "Let's setup the loop and run the simulations. This will take some time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "copyrighted-advance", + "metadata": {}, + "outputs": [], + "source": [ + "for exp in exp_list:\n", + " # Call run_model once for each experiment\n", + " edu_utils.run_model(exp[0], exp[1], exp[2], exp[3])" + ] + }, + { + "cell_type": "markdown", + "id": "exotic-frame", + "metadata": {}, + "source": [ + "## A look at the results\n", + "\n", + "When the simulations has finished running it is time to take a look at the results. Since we're keeping the results in `exp_list` we can again leverage the power of loops to save us writing some code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "located-cosmetic", + "metadata": {}, + "outputs": [], + "source": [ + "labels = ['Default', '-10%', '+10%']\n", + "fig, ax = plt.subplots(figsize=(12, 5))\n", + "for exp, label in zip(exp_list, labels):\n", + " # Get the data and plot it RESULTS are kept at the third spot, index 2\n", + " exp[2].surfMB.plot(ax=ax, label=label)\n", + "plt.legend(); " + ] + }, + { + "cell_type": "markdown", + "id": "excellent-allowance", + "metadata": {}, + "source": [ + "
\n", + "
\n", + " \n", + " Question: Can you explain why the surface mass balance looks like it does? Click me for an explanation!\n", + "

We are looking at three different surface mass balance series calculated with three different values for the albedo of fresh snow. The regular dips in the surface mass balance is the diurnal cycle, the glacier lose mass during the day while it remains constant during the night. When the albedo is increased the surface of the glacier will reflect more of the incoming solar radiation during the day. This leads to a decrease in the energy available for melt, thus reducing the melt and resulting in a less negative surface mass balance.\n", + "

\n", + "On the other hand, when the albedo is decreased the reflectivity of the glacier is reduced. This increases the energy available for melt which is why the surface mass balance is more negative.

\n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "amended-position", + "metadata": {}, + "source": [ + "*Here you can write your answer (double click to open the cell)* \n" + ] + }, + { + "cell_type": "markdown", + "id": "white-account", + "metadata": {}, + "source": [ + "## Another sensitivity study\n", + "\n", + "Try to come up with another variable to use for a sensitivity study. The loop for running the simulation is provided below, but you have to initialise the opt_dict(s), data and result datasets and the list containing your experiments." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "departmental-integration", + "metadata": {}, + "outputs": [], + "source": [ + "# Pick a variable to nudge\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "filled-french", + "metadata": {}, + "outputs": [], + "source": [ + "# Initialise the opt_dict(s) and datasets.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ordered-chassis", + "metadata": {}, + "outputs": [], + "source": [ + "# Pack it all into a list as above\n" + ] + }, + { + "cell_type": "markdown", + "id": "consistent-liberal", + "metadata": {}, + "source": [ + "
\n", + " Are you ready to run the simulation?\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "mechanical-titanium", + "metadata": {}, + "outputs": [], + "source": [ + "for exp in exp_list:\n", + " # Call run_model once for each experiment\n", + " edu_utils.run_model(exp[0], exp[1], exp[2], exp[3])" + ] + }, + { + "cell_type": "markdown", + "id": "electric-destination", + "metadata": {}, + "source": [ + "Try plotting the results of your study!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "through-rotation", + "metadata": {}, + "outputs": [], + "source": [ + "# Do some plotting here.\n" + ] + }, + { + "cell_type": "markdown", + "id": "thorough-material", + "metadata": {}, + "source": [ + "## Next steps\n", + "[Back to overview](welcome.ipynb)\n", + "\n", + "[Temperature bias experiments](temp_sensitivity.ipynb)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cosipy", + "language": "python", + "name": "cosipy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/temp_sensitivity.ipynb b/notebooks/temp_sensitivity.ipynb new file mode 100644 index 00000000..8a83641b --- /dev/null +++ b/notebooks/temp_sensitivity.ipynb @@ -0,0 +1,496 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "quantitative-appreciation", + "metadata": {}, + "source": [ + "# Input data and temperature sensitivity\n", + "\n", + "A continuation of the simple sensitivity study we did in \n", + "[Sensitivity studies using COSIPY](sensitivity_study.ipynb) is to investigate how changing the surface temperature, or any other input variable, affects the calculation of the surface mass balance. COSIPY reads the meteorological input data from a **netcdf** file during the run. This means that during each time step the model reads the meteorological data from the corresponding time step in the input file. The input file can be either \"1D\" for a point simulation or \"2D\" for a distributed simulation. The input file can be based on either observed or modeled data.\n", + "\n", + "In order to conduct the temperature sensitivity study (or any study) from scratch we need to do the following:\n", + "- Create the input file from observations.\n", + "- Load the input file into a xarray dataset.\n", + "- Create copies of this dataset and add the bias to the variable in question.\n", + "\n", + "With this done we can run the simulations. But before we do any of this, we can take a look at the input data for the Zhadang glacier, which is shipped with COSPIY for theses tutorials." + ] + }, + { + "cell_type": "markdown", + "id": "allied-filename", + "metadata": {}, + "source": [ + "**The standard imports**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "outer-calculation", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to change the cwd for the ipython session, otherwise COSIPY\n", + "# will look for things in the wrong places.\n", + "import os\n", + "import sys\n", + "# This is not really a good method, if cell is re run we end up in the\n", + "# wrong directory.\n", + "os.chdir('./../')\n", + "sys.path.append(os.getcwd())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "pursuant-exhibit", + "metadata": {}, + "outputs": [], + "source": [ + "from cosipy.utils import edu_utils\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ordinary-renewal", + "metadata": {}, + "outputs": [], + "source": [ + "# Have to tell matplotlib to plot inline\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "eligible-disorder", + "metadata": {}, + "source": [ + "The path to the input netcdf is stored in the `input_netcdf` variable. We can take a look at the options for the current value." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "concrete-wyoming", + "metadata": {}, + "outputs": [], + "source": [ + "# print options returns a pandas dataframe with the variable as the index,\n", + "# hence we can index it.\n", + "edu_utils.print_options().loc['input_netcdf']" + ] + }, + { + "cell_type": "markdown", + "id": "legislative-insulation", + "metadata": {}, + "source": [ + "We can then open the file with xarray:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fitted-start", + "metadata": {}, + "outputs": [], + "source": [ + "# Files are located in ./data/ relative to the root directory of cosipy. \n", + "# This is also our current working directory.\n", + "input_path = './data/' + 'input/' +\\\n", + " edu_utils.OPTIONS['input_netcdf']\n", + "with xr.open_dataset(input_path) as ds:\n", + " ds = ds.isel(time=slice(0, -1)).load()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "facial-yesterday", + "metadata": {}, + "outputs": [], + "source": [ + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "meaning-vertical", + "metadata": {}, + "source": [ + "As you can see, this file contain variables such as the 2-meter temperature (T2), the relative humidity at 2 meters (RH2) and the cloud cover (N) for one single point. The time coordinate spans 2009-01-01 00:00:00 and 2009-01-31 22:00:00 at a hourly resolution.\n", + "\n", + "\n", + "
\n", + "
\n", + " Question: Can you figure out what the variable G is and what unit it has? Click me for a hint!\n", + " Try pressing the document symbol to the right in the output above, this shows some extra information about the variable.\n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "instant-eleven", + "metadata": {}, + "source": [ + "We can quickly take a look at one of the variables by plotting it" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "turkish-dance", + "metadata": {}, + "outputs": [], + "source": [ + "ds.T2.plot(figsize=(10,5));" + ] + }, + { + "cell_type": "markdown", + "id": "alpine-debut", + "metadata": {}, + "source": [ + "## Creating the input files" + ] + }, + { + "cell_type": "markdown", + "id": "planned-february", + "metadata": {}, + "source": [ + "In this case COSIPY comes packaged with the processed data which can directly be used to drive a simulation. In a more realistic scenario however, you probably want drive COSIPY with your own data from another glacier than the Zhadang glacier. The `aws2cosipy` module is provided by COSIPY to aid the processing of .csv-files from weather stations or nwp-output into input files which can be used by COSIPY. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "freelance-selection", + "metadata": {}, + "outputs": [], + "source": [ + "# We import the create_1D_input function to process a single point\n", + "from utilities.aws2cosipy.aws2cosipy import create_1D_input" + ] + }, + { + "cell_type": "markdown", + "id": "corporate-middle", + "metadata": {}, + "source": [ + "The `create_1D_input` has five arguments: the csv file to process, the name of the resulting input file, the name of a static file describing the altitude, slope and aspect of the grid point, and the start and end date.\n", + "\n", + "In our case the .csv file is located in the same directory as the input file we looked at earlier. It is called `Zhadang_ERA5_2009_2018.csv`. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "irish-spare", + "metadata": {}, + "outputs": [], + "source": [ + "# Path to the file\n", + "file = './data/input/Zhadang/Zhadang_ERA5_2009_2018.csv'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "desirable-drilling", + "metadata": {}, + "outputs": [], + "source": [ + "# Static file\n", + "static_file = './data/static/Zhadang_static.nc'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "thick-fields", + "metadata": {}, + "outputs": [], + "source": [ + "# Define time start and end.\n", + "start_date = '20120701'\n", + "end_date = '20120731'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "hourly-laundry", + "metadata": {}, + "outputs": [], + "source": [ + "# Output file, we change this for our year.\n", + "output_name = './data/input/Zhadang/Zhadang_ERA5_2012_test.nc'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "raising-combination", + "metadata": {}, + "outputs": [], + "source": [ + "# Call the function, this takes some time.\n", + "create_1D_input(file, output_name, static_file, start_date, end_date)" + ] + }, + { + "cell_type": "markdown", + "id": "regular-phrase", + "metadata": {}, + "source": [ + "Open the file to see that it worked" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "opening-poland", + "metadata": {}, + "outputs": [], + "source": [ + "with xr.open_dataset(output_name) as ds:\n", + " ds = ds.isel(time=slice(0, -1)).load()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "regular-illness", + "metadata": {}, + "outputs": [], + "source": [ + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "noticed-lotus", + "metadata": {}, + "source": [ + "This created an input file with data for the month of July in 2012. The .csv-file contains data from January 2000 until early January 2018, so you are free to change the start and end data accordingly." + ] + }, + { + "cell_type": "markdown", + "id": "afraid-climate", + "metadata": {}, + "source": [ + "## Creating datasets for temperature bias experiments\n", + "\n", + "A temperature bias experiment is essentially a sensitivity study. It explores the effects of changing the temperature by $n$ degrees over the whole time period. For example, let's say we want to know how the glacier responds to a climate which is 2 °C warmer. This approach isolates the effect of temperature on the glacier and does not take other factors, like a changing hydrological cycle, into account.\n", + "\n", + "Setting this up with COSIPY follows a similar approach as described in the [Sensitivity studies with COSIPY](sensitivity_study.ipynb). However, instead of manipulating the the constants of the model with the `opt_dict` we are changing the dataset containing the input data before sending it to the `cosipy_core`.\n", + "\n", + "If we want to use our newly processed input file, we have to change the `input_netcdf` in the `opt_dict` along with the integration time (`time_start` & `time_end`) since this have to match the data. After this we can initiate the IOClass and datasets" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "quantitative-version", + "metadata": {}, + "outputs": [], + "source": [ + "opt_dict = {'input_netcdf': 'Zhadang/Zhadang_ERA5_2012_test.nc',\n", + " 'time_start': '2012-07-01T00:00',\n", + " 'time_end': '2012-07-31T00:00'}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "narrative-briefing", + "metadata": {}, + "outputs": [], + "source": [ + "# Default\n", + "IO_def, DATA_def, RESULTS_def = edu_utils.create_IO(opt_dict)\n", + "# Up\n", + "IO_up, DATA_up, RESULTS_up = edu_utils.create_IO(opt_dict)\n", + "# Down\n", + "IO_dn, DATA_dn, RESULTS_dn = edu_utils.create_IO(opt_dict)" + ] + }, + { + "cell_type": "markdown", + "id": "logical-elevation", + "metadata": {}, + "source": [ + "And then apply the bias to the input data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "suspended-charles", + "metadata": {}, + "outputs": [], + "source": [ + "# Bias is two degrees\n", + "bias = 2\n", + "DATA_dn['T2'] = DATA_dn['T2'] - bias\n", + "DATA_up['T2'] = DATA_up['T2'] + bias" + ] + }, + { + "cell_type": "markdown", + "id": "level-vietnam", + "metadata": {}, + "source": [ + "**Plot it to make sure it worked**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "baking-cement", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10,5))\n", + "DATA_def['T2'].plot(ax=ax, label='Default')\n", + "DATA_up['T2'].plot(ax=ax, label='+2')\n", + "DATA_dn['T2'].plot(ax=ax, label='-2')\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "continued-smith", + "metadata": {}, + "source": [ + "### Setting up the run\n", + "As in the previous sensitivity study we put everything in a list of lists. Note that we're not any of the options this time, only the input data and results dataset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cardiac-illinois", + "metadata": {}, + "outputs": [], + "source": [ + "# List of lists with our experiments\n", + "exp_list = [[DATA_def, IO_def, RESULTS_def],\n", + " [DATA_dn, IO_dn, RESULTS_dn],\n", + " [DATA_up, IO_up, RESULTS_up]\n", + " ]" + ] + }, + { + "cell_type": "markdown", + "id": "amended-interpretation", + "metadata": {}, + "source": [ + "And lets runs the experiments" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "pleasant-billion", + "metadata": {}, + "outputs": [], + "source": [ + "for exp in exp_list:\n", + " # Call run_model once for each experiment\n", + " edu_utils.run_model(DATA=exp[0], IO=exp[1], RESULT=exp[2])" + ] + }, + { + "cell_type": "markdown", + "id": "copyrighted-bikini", + "metadata": {}, + "source": [ + "We select a few days in the middle of the month to plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "radical-liver", + "metadata": {}, + "outputs": [], + "source": [ + "labels = ['Default', '-2$\\degree$C', '+2$\\degree$C']\n", + "fig, ax = plt.subplots(figsize=(10, 5))\n", + "for exp, label in zip(exp_list, labels):\n", + " # Get the data and plot it RESULTS are kept at the third spot, index 2\n", + " exp[2].sel(time=slice('2012-07-07', '2012-07-12')).surfMB.plot(ax=ax,\n", + " label=label)\n", + "plt.legend(); " + ] + }, + { + "cell_type": "markdown", + "id": "dangerous-salem", + "metadata": {}, + "source": [ + "
\n", + "
\n", + " \n", + " Question: What is the difference between the experiments in the total mass lost during the period plotted above? Click me for an explanation\n", + " \n", + " Try selecting the data from the experiment list and then calculate the sum with .sum(). Alternatively you can also copy the cell above and add .cumsum() before .plot(). This will plot the cumulative sum over the period.\n", + "
\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "sorted-future", + "metadata": {}, + "outputs": [], + "source": [ + "# Here you can write some code. Click to open.\n" + ] + }, + { + "cell_type": "markdown", + "id": "harmful-remedy", + "metadata": {}, + "source": [ + "## Next steps\n", + "[Back to overview](welcome.ipynb)\n", + "\n", + "[Distributed simulations](distributed_run.ipynb)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/welcome.ipynb b/notebooks/welcome.ipynb new file mode 100644 index 00000000..83d1e82f --- /dev/null +++ b/notebooks/welcome.ipynb @@ -0,0 +1,55 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![](https://cryo-tools.org/wp-content/uploads/2019/11/COSIPY-logo-2500px.png)\n", + "\n", + "# Welcome to the COSIPY educational notebooks!\n", + "\n", + "COSIPY - the **Co**upled **S**nowpack and **I**ce surface energy and mass balance model in **Py**thon - is an open-source model for simulating distributed snow and glacier mass changes.\n", + "\n", + "In the following notebooks we will go through some basic concepts of COSIPY, how it works and what you can do with it. \n", + "\n", + "If jupyter notebooks are a new concept to you, we recommend you to take a look at [this great introduction](https://edu.oggm.org/en/latest/notebooks_howto.html) from the [OGGM-team](https://oggm.org/).\n", + "\n", + "These notebooks can be rendered as a static webpage, or if you want to be able to play around with the model, open them on [MyBinder](https://mybinder.org/) by clicking the button below\n", + "\n", + "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/Holmgren825/cosipy/edu_revamp?filepath=notebooks%2Fwelcome.ipynb)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A list of the tutorials\n", + "- [First steps with COSIPY](first_steps.ipynb) Explains how to setup and run COSIPY on one grid point on a glaicer.\n", + "- [Sensitivity studies using COSIPY](sensitivity_study.ipynb) Show you how to do a simple sensitivity study using COSIPY.\n", + "- [Temperature bias experiments](temp_sensitivity.ipynb) Extends on the sensitivity study notebook and explains how to set up experiments to test the temperature sensitivity.\n", + "- [A distributed run](distributed_run.ipynb) Goes through how to set up a distributed run and briefly discusses how this work." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/requirements.txt b/requirements.txt index 64f737eb..a1adf7b2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,8 +3,9 @@ pandas numba scipy xarray<=0.18.2 -distributed<=2.5.2 +distributed dask_jobqueue netcdf4 coveralls metpy +matplotlib