{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectral filtering (spectral round-trip at a lower truncation) for global fields\n",
    "\n",
    "and vizualisation of differences as variance spectrum plot. \n",
    "Except for the truncation part, the approach proposed here using `H2DField.spectrum()` could also be used for limited-area fields. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-13T11:34:00.159904Z",
     "iopub.status.busy": "2025-01-13T11:34:00.158570Z",
     "iopub.status.idle": "2025-01-13T11:34:00.688149Z",
     "shell.execute_reply": "2025-01-13T11:34:00.687560Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# for figures in notebook\n",
    "\n",
    "# import & initialize epygram\n",
    "import epygram\n",
    "\n",
    "epygram.init_env()\n",
    "\n",
    "import os\n",
    "\n",
    "INPUTS_DIR = os.path.join(\"..\", \"inputs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-13T11:34:00.690407Z",
     "iopub.status.busy": "2025-01-13T11:34:00.690202Z",
     "iopub.status.idle": "2025-01-13T11:34:02.897233Z",
     "shell.execute_reply": "2025-01-13T11:34:02.896608Z"
    }
   },
   "outputs": [],
   "source": [
    "res = epygram.open(\n",
    "    os.path.join(INPUTS_DIR, \"bgstderr.arpege-unbal-ens.tl224+0003:00.grib\"), \"r\"\n",
    ")\n",
    "field = res.readfield(\"shortName=vo,level=1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some specific cases, the spectral geometry of the field cannot be accessed from the field itself. For GRIB files for instance, it may have to be deduced from the grid-point geometry, in the case of global Gauss geometries. This can be done using the `epygram.spectra.get_spectral_geometry()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deduce spectral geometry of GRIB field from grid-point geometry\n",
    "spgeom = epygram.spectra.get_spectral_geometry(field, res)\n",
    "if spgeom is None:\n",
    "    spgeom = epygram.spectra.make_spectral_geometry(field.geometry)\n",
    "\n",
    "# initial field spectrum with specified name\n",
    "sp01 = field.spectrum(\"Vorticity @Lev1\", spectral_geometry=spgeom)\n",
    "field.sp2gp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-13T11:34:03.089076Z",
     "iopub.status.busy": "2025-01-13T11:34:03.088811Z",
     "iopub.status.idle": "2025-01-13T11:34:03.213998Z",
     "shell.execute_reply": "2025-01-13T11:34:03.213475Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set new spectral space\n",
    "print(spgeom.truncation)\n",
    "new_truncation = spgeom.truncation.copy()\n",
    "new_truncation[\"max\"] = spgeom.truncation[\"max\"] // 2\n",
    "print(new_truncation)\n",
    "trunc_spgeom = epygram.geometries.SpectralGeometry(\"legendre\", new_truncation)\n",
    "# Spectral roundtrip\n",
    "field.gp2sp(trunc_spgeom)\n",
    "field.sp2gp()\n",
    "# And compute filtered spectrum\n",
    "sp01f = field.spectrum(sp01.name + \" (filtered)\", spectral_geometry=spgeom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-01-13T11:34:03.215956Z",
     "iopub.status.busy": "2025-01-13T11:34:03.215736Z",
     "iopub.status.idle": "2025-01-13T11:34:03.678853Z",
     "shell.execute_reply": "2025-01-13T11:34:03.678243Z"
    }
   },
   "outputs": [],
   "source": [
    "# and plot\n",
    "fig, ax = epygram.spectra.plotspectra([sp01, sp01f])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}