{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A Pseudoproxy Experiment with GraphEM" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Expected time to run through: 10 mins**\n", "\n", "This tutorial demonstrates how to get a reconstruction using GraphEM, leveraging a simple pseudoproxy dataset generated from a subset of iCESM gridded dataset.\n", "The pseudoproxy are generated based on the original iCESM simulated `tas` plus white noise with `SNR=1`, using below code block:\n", "\n", "```python\n", "# gen pseudoproxy\n", "\n", "SNR = 1\n", "\n", "pp = np.copy(tas_sub)\n", "nt, nlat, nlon = np.shape(pp)\n", "\n", "k = 0\n", "for i in range(nlat):\n", " for j in range(nlon):\n", " tas_std = np.std(tas_sub[:,i,j])\n", " noise_std = tas_std / SNR\n", " \n", " np.random.seed(k)\n", " pp[:,i,j] += np.random.normal(0, noise_std, size=nt)\n", " \n", " k += 1\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test data preparation\n", "\n", "To go through this tutorial, please prepare test data following the steps:\n", "\n", "Download the test case named “PPE” with this [link](https://drive.google.com/drive/folders/1QpXqFgCBA4QMb7vqGIWxptkhWTJp0e8c?usp=sharing).\n", "Create a directory named “testcases” in the same directory where this notebook sits.\n", "Put the unzipped direcotry “PPE” into “testcases”.\n", "\n", "Below, we first load some useful packages, including our `GraphEM`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import LMRt\n", "import GraphEM\n", "import os\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reconstruction" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "job = GraphEM.ReconJob()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.load_configs() >>> loading reconstruction configurations from: ./testcases/PPE/configs.yml\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_configs() >>> job.configs created\u001b[0m\n", "\u001b[1m\u001b[36mGraphEM: job.load_configs() >>> job.configs[\"job_dirpath\"] = /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PPE/recon\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_configs() >>> /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PPE/recon created\u001b[0m\n", "{'anom_period': [1951, 1980],\n", " 'calib_period': [1750, 1849],\n", " 'job_dirpath': '/Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PPE/recon',\n", " 'job_id': 'LMRt_quickstart',\n", " 'obs_path': {'tas': './data/obs/iCESM_subset.nc'},\n", " 'obs_varname': {'lat': 'lat', 'lon': 'lon', 'tas': 'tas'},\n", " 'proxydb_path': './data/proxy/pseudoproxy_dataset.pkl',\n", " 'ptype_list': ['coral.d18O'],\n", " 'recon_period': [850, 1849]}\n" ] } ], "source": [ "job.load_configs('./testcases/PPE/configs.yml', verbose=True)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.load_proxydb() >>> job.configs[\"proxydb_path\"] = /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PPE/data/proxy/pseudoproxy_dataset.pkl\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_proxydb() >>> 100 records loaded\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.load_proxydb() >>> job.proxydb created\u001b[0m\n" ] } ], "source": [ "job.load_proxydb(verbose=True)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.filter_proxydb() >>> filtering proxy records according to: ['coral.d18O']\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.filter_proxydb() >>> 100 records remaining\u001b[0m\n" ] } ], "source": [ "job.filter_proxydb(verbose=True)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[36mGraphEM: job.load_obs() >>> loading instrumental observation fields from: {'tas': '/Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PPE/data/obs/iCESM_subset.nc'}\u001b[0m\n", "Time axis not overlap with the reference period [1951, 1980]; use its own time period as reference [850.00, 1849.00].\n", "\u001b[1m\u001b[32mGraphEM: job.load_obs() >>> job.obs created\u001b[0m\n" ] } ], "source": [ "job.load_obs(verbose=True)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset Overview\n", "-----------------------\n", "\n", " Name: tas\n", " Source: /Users/fzhu/Github/GraphEM/docsrc/tutorial/testcases/PPE/data/obs/iCESM_subset.nc\n", " Shape: time:1000, lat:10, lon:10\n" ] } ], "source": [ "print(job.obs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the loaded iCESM simulation and the pseudoproxy dataset are already annualized, we can skip the `.seasonalize()` steps and run `.prep_data()` directly." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.recon_time created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.calib_time created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.calib_idx created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.temp created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.df_proxy created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.proxy created\u001b[0m\n", "\u001b[1m\u001b[32mGraphEM: job.prep_data() >>> job.lonlat created\u001b[0m\n" ] } ], "source": [ "job.prep_data(verbose=True)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | pp_000 | \n", "pp_001 | \n", "pp_002 | \n", "pp_003 | \n", "pp_004 | \n", "pp_005 | \n", "pp_006 | \n", "pp_007 | \n", "pp_008 | \n", "pp_009 | \n", "... | \n", "pp_090 | \n", "pp_091 | \n", "pp_092 | \n", "pp_093 | \n", "pp_094 | \n", "pp_095 | \n", "pp_096 | \n", "pp_097 | \n", "pp_098 | \n", "pp_099 | \n", "
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
850.0 | \n", "300.751783 | \n", "300.699385 | \n", "299.718651 | \n", "300.717691 | \n", "299.793416 | \n", "299.939352 | \n", "299.472350 | \n", "300.551011 | \n", "299.553801 | \n", "299.413664 | \n", "... | \n", "299.249156 | \n", "299.027787 | \n", "299.182322 | \n", "299.526141 | \n", "299.161496 | \n", "298.628304 | \n", "298.802706 | \n", "298.107757 | \n", "298.818966 | \n", "298.544758 | \n", "
851.0 | \n", "300.496749 | \n", "300.009875 | \n", "300.193606 | \n", "300.365040 | \n", "300.347096 | \n", "299.853183 | \n", "300.380611 | \n", "299.650771 | \n", "300.502845 | \n", "299.614965 | \n", "... | \n", "299.103565 | \n", "298.778987 | \n", "299.141563 | \n", "299.335989 | \n", "299.410467 | \n", "298.670001 | \n", "298.838717 | \n", "298.689713 | \n", "299.370283 | \n", "299.494514 | \n", "
852.0 | \n", "300.902166 | \n", "300.334244 | \n", "299.602851 | \n", "300.671260 | \n", "300.102458 | \n", "301.886292 | \n", "300.641600 | \n", "300.462630 | \n", "299.210613 | \n", "299.595454 | \n", "... | \n", "299.289380 | \n", "299.249218 | \n", "299.134434 | \n", "299.942948 | \n", "298.988984 | \n", "299.277950 | \n", "298.905532 | \n", "298.840735 | \n", "299.414093 | \n", "299.140887 | \n", "
853.0 | \n", "301.729269 | \n", "300.347924 | \n", "301.572955 | \n", "299.850196 | \n", "301.104819 | \n", "300.567805 | \n", "300.145062 | \n", "300.817268 | \n", "299.709402 | \n", "300.448917 | \n", "... | \n", "299.573146 | \n", "299.649931 | \n", "298.838391 | \n", "299.710779 | \n", "300.040319 | \n", "299.609617 | \n", "299.473822 | \n", "299.656169 | \n", "299.019936 | \n", "299.719184 | \n", "
854.0 | \n", "301.378359 | \n", "300.997108 | \n", "299.769165 | \n", "300.429641 | \n", "300.280460 | \n", "300.495473 | \n", "298.998074 | \n", "299.858000 | \n", "298.895901 | \n", "299.967372 | \n", "... | \n", "300.078972 | \n", "300.040279 | \n", "299.551374 | \n", "300.088576 | \n", "298.790187 | \n", "299.790875 | \n", "299.441056 | \n", "299.184892 | \n", "299.035147 | \n", "299.345002 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
1845.0 | \n", "300.370235 | \n", "300.113197 | \n", "300.144971 | \n", "299.726555 | \n", "299.921271 | \n", "300.382517 | \n", "299.913085 | \n", "299.789860 | \n", "299.300107 | \n", "299.261318 | \n", "... | \n", "299.482236 | \n", "299.063248 | \n", "299.152310 | \n", "298.784250 | \n", "298.994037 | \n", "298.300974 | \n", "299.091353 | \n", "298.409069 | \n", "298.584221 | \n", "298.305443 | \n", "
1846.0 | \n", "299.985712 | \n", "298.994364 | \n", "300.207191 | \n", "299.883894 | \n", "300.684450 | \n", "300.181971 | \n", "300.745885 | \n", "299.919197 | \n", "298.775608 | \n", "299.324020 | \n", "... | \n", "299.544706 | \n", "299.056843 | \n", "299.350263 | \n", "299.165946 | \n", "299.459059 | \n", "299.227588 | \n", "299.156746 | \n", "299.367724 | \n", "299.111197 | \n", "298.400678 | \n", "
1847.0 | \n", "300.292753 | \n", "300.177510 | \n", "300.033955 | \n", "300.236235 | \n", "299.559872 | \n", "299.877648 | \n", "299.141349 | \n", "299.497879 | \n", "299.952182 | \n", "299.577197 | \n", "... | \n", "299.724081 | \n", "300.424495 | \n", "299.179294 | \n", "299.375003 | \n", "299.475237 | \n", "299.947260 | \n", "300.159204 | \n", "299.223168 | \n", "299.384652 | \n", "299.483399 | \n", "
1848.0 | \n", "299.823016 | \n", "300.400975 | \n", "300.223044 | \n", "300.913630 | \n", "300.261049 | \n", "300.159285 | \n", "299.311807 | \n", "299.449942 | \n", "299.660750 | \n", "298.979497 | \n", "... | \n", "299.829125 | \n", "299.592597 | \n", "300.013318 | \n", "299.155958 | \n", "299.107011 | \n", "299.086289 | \n", "299.878042 | \n", "299.307766 | \n", "299.593688 | \n", "298.706442 | \n", "
1849.0 | \n", "300.161352 | \n", "300.159181 | \n", "300.246488 | \n", "299.820504 | \n", "299.517483 | \n", "299.951131 | \n", "299.332212 | \n", "299.316922 | \n", "299.630661 | \n", "299.881939 | \n", "... | \n", "299.587589 | \n", "298.878817 | \n", "298.997188 | \n", "298.895755 | \n", "299.134630 | \n", "299.284257 | \n", "299.369434 | \n", "299.400078 | \n", "299.026855 | \n", "299.096824 | \n", "
1000 rows × 100 columns
\n", "