Step1a_Generate_RAS_fixed_resolution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import dask\n",
    "import dask.array as da\n",
    "from dask.diagnostics import ProgressBar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearsim import time_domain_ras"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate a dataset of input spectra "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_spec = xr.Dataset()\n",
    "ds_spec['Gamma'] = xr.DataArray([3.3],dims='gamma')\n",
    "ds_spec['Tp'] = xr.DataArray([12.],dims='tp')\n",
    "ds_spec['Duration'] = xr.DataArray(da.array([2,5,10,20,40,60])*60,dims='duration')\n",
    "n_seed = 1000\n",
    "ds_spec['Seed'] = xr.DataArray(da.arange(n_seed,chunks=(200,)).astype('int'),dims='seed')\n",
    "ds_spec['this_seed'] = xr.DataArray(da.random.choice(int(1E9), size=(len(ds_spec['Tp']), len(ds_spec['Gamma']), len(ds_spec['Duration']), n_seed)).astype('int'), dims=('tp','gamma','duration','seed'))\n",
    "ds_spec = ds_spec.chunk({'tp':1,'gamma':-1,'duration':1,'seed':200})\n",
    "ds_spec['hs'] = 1.0\n",
    "ds_spec['dt'] = 0.5\n",
    "ds_spec['fft_min_duration']=60*60\n",
    "\n",
    "def set_coords(ds):\n",
    "    ds = ds.rename({'Duration':'duration','Gamma':'gamma','Seed':'seed','Tp':'tp'})\n",
    "    ds = ds.assign_coords({'duration':ds['duration'],'gamma':ds['gamma'],'seed':range(len(ds['seed'])),'tp':ds['tp']})\n",
    "    ds = ds.set_coords(['hs','tp','dt','fft_min_duration'])\n",
    "    return ds\n",
    "\n",
    "ds_spec = set_coords(ds_spec)\n",
    "\n",
    "ds_spec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Names of output variables returned from time_domain_ras function\n",
    "output_variables = 'Tz, Tm01, Hm0, Hs, H13, H13_unbiased, Hmax, HmaxT, Cmax, CmaxT2, r_spectra, r_sample, r_unbiased, k3, k4, seed'\n",
    "outputs = xr.apply_ufunc(time_domain_ras,\n",
    "                        ds_spec['tp'],\n",
    "                        ds_spec['hs'],\n",
    "                        ds_spec['gamma'],\n",
    "                        ds_spec['duration'],\n",
    "                        ds_spec['dt'],\n",
    "                        ds_spec['this_seed'],\n",
    "                        ds_spec['fft_min_duration'],\n",
    "                        input_core_dims=[[],[],[],[],[],[],[]],\n",
    "                        output_core_dims=[[],]*16,\n",
    "                        vectorize=True,\n",
    "                        dask='parallelized',\n",
    "                        output_dtypes=['float',]*16\n",
    "                        )\n",
    "                        \n",
    "# Assign the outputs to the xr.Dataset\n",
    "for v, o in zip(output_variables.split(',')[:-1],outputs[:-1]): # ignore the returned seed - same as input\n",
    "    ds_spec[v.strip()] = o\n",
    "ds_spec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use multiple processes to calculate timeseries\n",
    "with dask.config.set(scheduler='processes'):\n",
    "    with ProgressBar():\n",
    "        ds_spec = ds_spec.compute()\n",
    "ds_spec.to_netcdf(f'data/timeseries_stats/RAS_n{n_seed}_2Hz_fixedres.nc',mode='w')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyse saved results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_spec = xr.open_dataset(f'data/timeseries_stats/RAS_n{n_seed}_2Hz_fixedres.nc')\n",
    "ds_spec['HmHs'] = ds_spec['Hmax']/ds_spec['Hs']\n",
    "ds_std = ds_spec.std(dim='seed')\n",
    "ds_mean = ds_spec.mean(dim='seed')\n",
    "ds_max = ds_spec.max(dim='seed')\n",
    "ds_CV = ds_std/ds_mean\n",
    "ds_CV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.xarray\n",
    "import holoviews as hv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layout = []\n",
    "for v in ds_mean.data_vars:\n",
    "    layout.append(ds_mean[v].hvplot(x='duration',cmap='viridis',aspect=1,title=v))\n",
    "layout = hv.Layout(layout).cols(3)\n",
    "\n",
    "hv.save(layout,filename=\"figures/timeseries_stats/RAS_Mean_fixedres.html\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layout = []\n",
    "for v in ds_CV.data_vars:\n",
    "    layout.append(ds_CV[v].hvplot(x='duration',cmap='viridis',aspect=1,title=v))\n",
    "layout = hv.Layout(layout).cols(3)\n",
    "\n",
    "hv.save(layout,filename=\"figures/timeseries_stats/RAS_CV_fixedres.html\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layout = []\n",
    "for v in ds_max.data_vars:\n",
    "    layout.append(ds_max[v].hvplot(x='duration',cmap='viridis',aspect=1,title=v))\n",
    "layout = hv.Layout(layout).cols(3)\n",
    "\n",
    "hv.save(layout,filename=\"figures/timeseries_stats/RAS_Max_fixedres.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take a look at the largest wave event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_stacked = ds_spec.stack(sample=['tp','duration','gamma','seed'])\n",
    "im=ds_stacked.Hmax.argmax()\n",
    "largest_wave = ds_stacked.isel(sample=im)\n",
    "\n",
    "display(largest_wave)\n",
    "\n",
    "tp, duration, gamma,_ = np.array(largest_wave.sample.values.item())\n",
    "seed = int(largest_wave.this_seed)\n",
    "dt=float(ds_stacked['dt'])\n",
    "hs=float(ds_stacked['hs'])\n",
    "fft_duration=ds_stacked['fft_min_duration']\n",
    "ts=time_domain_ras(tp,hs,gamma,duration,dt=dt,seed=seed,return_ts=True,fft_equiv_duration=fft_duration)\n",
    "da_ts = xr.DataArray(ts,dims='time',coords={'time':np.arange(0,duration,dt)},name='timeseries')\n",
    "\n",
    "# Double check that that the time series is generates correctly\n",
    "from linearsim import wave_stats\n",
    "print(f'ts statistics: {wave_stats(ts,1/dt)}')\n",
    "\n",
    "da_ts.hvplot(x='time')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('spectem')",
   "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.9.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c44064dc6b1eba8586ad4f01501b43f80852b132e8c1f3c24e235fa6d61993b5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}