{ "cells": [ { "cell_type": "markdown", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "# Calibrating GPs using ABC" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings('ignore') # Ignore all the iris warnings..." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import cis\n", "import iris\n", "\n", "from utils import get_aeronet_data, get_bc_ppe_data\n", "\n", "from esem.utils import validation_plot, plot_parameter_space, get_random_params, ensemble_collocate\n", "from esem import gp_model\n", "from esem.abc_sampler import ABCSampler, constrain\n", "\n", "import iris.quickplot as qplt\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in the parameters and observables" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ungridded data: Absorption_AOD440nm / (1) \n", " Shape = (10098,)\n", "\n", " Total number of points = 10098\n", " Number of non-masked points = 10098\n", " Long name = Absorption_AOD440nm\n", " Standard name = None\n", " Units = 1\n", " Missing value = -999.0\n", " Range = (5.1e-05, 0.47236)\n", " History = \n", " Coordinates: \n", " longitude\n", " Long name = \n", " Standard name = longitude\n", " Units = degrees_east\n", " Missing value = None\n", " Range = (-155.576755, 141.3407)\n", " History = \n", " latitude\n", " Long name = \n", " Standard name = latitude\n", " Units = degrees_north\n", " Missing value = None\n", " Range = (-35.495807, 79.990278)\n", " History = \n", " time\n", " Long name = \n", " Standard name = time\n", " Units = days since 1600-01-01 00:00:00\n", " Missing value = None\n", " Range = (cftime.DatetimeGregorian(2017, 1, 1, 12, 0, 0, 0, has_year_zero=False), cftime.DatetimeGregorian(2018, 1, 2, 12, 0, 0, 0, has_year_zero=False))\n", " History = \n", "\n" ] } ], "source": [ "aaod = get_aeronet_data()\n", "print(aaod)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Read in the PPE parameters, AAOD and DRE\n", "ppe_params, ppe_aaod, ppe_dre = get_bc_ppe_data(dre=True)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:root:Creating guessed bounds as none exist in file\n", "WARNING:root:Creating guessed bounds as none exist in file\n", "WARNING:root:Creating guessed bounds as none exist in file\n", "WARNING:root:Creating guessed bounds as none exist in file\n" ] } ], "source": [ "# Take the annual mean of the DRE\n", "ppe_dre, = ppe_dre.collapsed('time', iris.analysis.MEAN)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Collocate the model on to the observations" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "col_ppe_aaod = ensemble_collocate(ppe_aaod, aaod)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "n_test = 8\n", "\n", "X_test, X_train = ppe_params[:n_test], ppe_params[n_test:]\n", "Y_test, Y_train = col_ppe_aaod[:n_test], col_ppe_aaod[n_test:]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
Absorption Optical Thickness - Total 550Nm (1) | \n", "job | \n", "obs | \n", "
---|---|---|
Shape | \n", "31 | \n", "10098 | \n", "
Dimension coordinates | \n", "\n", " | \n", " |
\tjob | \n", "x | \n", "- | \n", "
\tobs | \n", "- | \n", "x | \n", "
Auxiliary coordinates | \n", "\n", " | \n", " |
\tlatitude | \n", "- | \n", "x | \n", "
\tlongitude | \n", "- | \n", "x | \n", "
\ttime | \n", "- | \n", "x | \n", "