paper_videos.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af25a731",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%javascript #enables smooth scrolling in jupyter notebook\n",
    "Jupyter.keyboard_manager.command_shortcuts.remove_shortcut('up');\n",
    "Jupyter.keyboard_manager.command_shortcuts.remove_shortcut('down');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50632853",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "display(HTML(\"<style>.rendered_html { font-size: 20px; }</style>\"))\n",
    "\n",
    "import sys\n",
    "#sys.path.append('/home/mchoblet/paleoda')\n",
    "\n",
    "from global_land_mask import globe\n",
    "import matplotlib.gridspec as gridspec\n",
    "import xskillscore as xss\n",
    "from os import listdir\n",
    "import os\n",
    "\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import string\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature  \n",
    "from matplotlib import cm as cm\n",
    "import numpy as np\n",
    "import nc_time_axis #for plotting time axis correctly with .plot() (had to update that + xarray to newest version)\n",
    "import scipy.stats as stats\n",
    "\n",
    "import matplotlib as mpl\n",
    "import pyleoclim as pyleo\n",
    "mpl.rc('hatch', color='darkgray', linewidth=1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcd1d701",
   "metadata": {},
   "source": [
    "# Converting reconstruction plots into videos\n",
    "\n",
    "**Procedure adopted here:**\n",
    "* creation of png snapshots which are saved in numerated order. So basically its all about creating many plots.\n",
    "* The pngs concatenated into an animation using _ffmpeg_ command line tool (this is also where repetion rate is set (see description below)\n",
    "* I choose a coarse temporal updating here (25/50 years), though the reconstruction data is annually (but I don't want too many numbers wiggling around)\n",
    "\n",
    "**Things to take into account:**\n",
    "* It is crucial to set the position of elements (axis limits, axis labels, legends...) explicitely such that they don't change between the different pngs\n",
    "* There are also python animation tools to create the video directly, however I struggled a bit and choose the save way why pngs\n",
    "* Setting the time axis labels to specific years is a bit annoying (but the result is acceptable)\n",
    "\n",
    "**Videos produced here**:\n",
    "* Climate field anomalies during the Common Era\n",
    "* Precipitation field + Monsoon Index + Neukom 2018 GMST and GMST core monsoon temperature \n",
    "* Precipitation field + Monsoon Index + Speleothem anomalies + Neukom 2018 GMST and GMST core monsoon temperature\n",
    "\n",
    "Note, that while only speleothem anomalies are displayed, also other data influence the reconstruction. That data uses other units, so putting everything into a single plot would be messy. You can also use the code here to produce an animation of the reconstruction which only relies on speleothem data.\n",
    "\n",
    "**Data used for comparison**\n",
    "* Neukom, Raphael; Steiger, Nathan; Gómez-Navarro, Juan José; Werner, Johannes; Wang, Jianghao (2019). No evidence for globally coherent warm and cold periods over the pre-industrial Common Era. figshare. Collection. https://doi.org/10.6084/m9.figshare.c.4498373 \n",
    "\n",
    "This is kind of the most widely accepted temperature reconstruction of the Common Era (It's a multi-method mean). I don't use the data in the manuscript elsewhere (there I focus on the LMRv2.1. and PHYDA reconstruction), because this dataset is only for temperature. The main purpose of the paper is also not the comparison to other datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b96e12a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def anomaly(ds,start='1800',end='1850'):\n",
    "    \"\"\"\n",
    "    start='1850'\n",
    "    end='1900'\n",
    "    \"\"\"\n",
    "    ds=ds-ds.sel(time=slice(start,end)).mean('time')\n",
    "    return ds\n",
    "\n",
    "def quick_filt(ts,cutoff=1/30):\n",
    "    from scipy.signal import butter, lfilter, filtfilt\n",
    "    def butter_lowpass(cutoff, fs, order=6, kf=1.2):\n",
    "        # kf:  scaling factor for the lowpass frequency; 1 = Nyquist, 1.2 =\n",
    "        #' 1.2xNyquist is a tradeoff between reducing variance loss and keeping\n",
    "        #' aliasing small\n",
    "        #fs is basic timestep (min_ts)\n",
    "        #nyquist frequency\n",
    "        nyq = 0.5 * fs \n",
    "        normal_cutoff = cutoff / nyq * kf\n",
    "        b, a = butter(order, normal_cutoff, btype='low', analog=False)\n",
    "        return b, a\n",
    "    def butter_lowpass_filter(data, cutoff, fs, order=6):\n",
    "        #filtfilt does not introduce a time-lag in comparison to butterworth \n",
    "        b, a = butter_lowpass(cutoff, fs, order=order)\n",
    "        y = filtfilt(b, a, data)\n",
    "        #y = lfilter(b, a, data)\n",
    "        return y\n",
    "    return butter_lowpass_filter(ts,cutoff,fs=1,order=4)\n",
    "\n",
    "monsoon_region=[[-17.5,-5],[287.5,312.5]]\n",
    "\n",
    "def regionalmean(field,lats=monsoon_region[0],lons=monsoon_region[1],name=None):\n",
    "    \"\"\"\n",
    "    Function that calculates the regional mean of a climate field. (Eg nino index if you have SST data).\n",
    "    Does latitudinal weight averageing. Also computes correct average for cross zero meridan regions (e.g 350,10)\n",
    "\n",
    "    Input:\n",
    "        - Climate Field (time,lat,lon): GMT computed over lat,lon for given restructions\n",
    "        - lat =[latSouth,latNorth], lon =[latWest,lonWest] Limits of the Region of interest.\n",
    "        Given in -90 to 90 and 0 to 360 degrees.\n",
    "        - variable name\n",
    "\n",
    "    Maybe in the future: Non-rectangular regions (there I would need to explicitely specify all the gridpoints I guess)\n",
    "    \"\"\"\n",
    "    lat=field.lat\n",
    "    lon=field.lon\n",
    "    \n",
    "    #selection of longitudes\n",
    "    if lons[0]<lons[1]:\n",
    "        sel_lat = lat.where( (lat >= lats[0] ) & (lat <= lats[1]), drop=True)\n",
    "        sel_lon = lon.where( (lon >= lons[0] ) & (lon <= lons[1]), drop=True)\n",
    "    #cross zero-meridian region option\n",
    "    else:\n",
    "        sel_lat = lat.where( (lat >= lats[0] ) & (lat <= lats[1]), drop=True)\n",
    "        sel_lon_1 = lon.where( (lon >= lons[0] ), drop=True)\n",
    "        sel_lon_2 = lon.where((lon <= lons[1]), drop=True)\n",
    "        sel_lon=xr.concat([sel_lon_1,sel_lon_2],dim='lon')\n",
    "    \n",
    "    wgt=np.cos(np.deg2rad(lat))\n",
    "    field_m=field.sel(lat=sel_lat,lon=sel_lon).weighted(wgt).mean(('lat','lon'))\n",
    "    if name!=None:\n",
    "        field_m=field_m.rename(('regm_'+name))\n",
    "    return field_m\n",
    "\n",
    "def globalmean(field,name=''):\n",
    "    \"\"\"\n",
    "    Function that calculates the global mean of a climate field \n",
    "    (DataArray from one variable, not all variables)\n",
    "    Does latitude weighting.\n",
    "\n",
    "    Input:\n",
    "        - Climate Field (time,lat,lon): GMT computed over lat,lon\n",
    "        - variable name as string, (needed for naming)\n",
    "    Output:\n",
    "        - Global mean as DataArray with name \"globalmean_+<name>\"\n",
    "        \n",
    "    Other possible indices\n",
    "    - El -Nino index: Sea Surface Temperature(5°N - 5°S), 120° - 170° West (190-240° in 0-360) mode\n",
    "    (https://en.wikipedia.org/wiki/Multivariate_ENSO_index) (I currently don't have SST)\n",
    "    - AMOC index: usually defined as the stream function for the zonally integrated meridional volume transport in the ocean\n",
    "    (at some specific latitude). Steiger 2016: (defined here as the maximum value of the overturn-\n",
    "    ing stream function in the North Atlantic between 25 and\n",
    "    70◦ N and between depths of 500 and 2000 m\n",
    "    - Convert -180,180 to 0,360 -> lon % 360\n",
    "    \"\"\"\n",
    "    lat=field.lat\n",
    "    wgt=np.cos(np.deg2rad(lat))\n",
    "    field_m=field.weighted(wgt).mean(('lat','lon'))\n",
    "    try:\n",
    "        field_m=field_m.rename(('gm_'+name))\n",
    "    except:\n",
    "        pass\n",
    "    return field_m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83d39441",
   "metadata": {},
   "source": [
    "# Load Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5955f67",
   "metadata": {},
   "outputs": [],
   "source": [
    "### load the reconstruction data and resample it to Decadal Resolution\n",
    "folder='/home/mchoblet/paleoda/results/experiments/paper/revision/normal_err/annual/proxy_frac/all_proxies/'\n",
    "variables=['tsurf_mean','prec_mean','spei_mean','d18O_mean']\n",
    "files=listdir(folder)\n",
    "file_l=[xr.open_dataset(folder+f,use_cftime=True)[variables] for f in files]\n",
    "combined_dataset = xr.concat(file_l, dim='dataset')\n",
    "recon=combined_dataset.mean('dataset')\n",
    "recon=recon-recon.sel(time=slice('0850','1850')).mean('time') #demean!\n",
    "recon['spei_mean']=(recon['spei_mean']-recon['spei_mean'].sel(time=slice('0850','1850')).mean('time'))/recon['spei_mean'].sel(time=slice('0850','1850')).std('time')\n",
    "\n",
    "#limit spei to values on earth\n",
    "data=recon['spei_mean']\n",
    "data['lon']=(data.lon+180)%360-180\n",
    "data=data.sortby('lon')\n",
    "lat=data['lat']\n",
    "lon=data['lon']\n",
    "lat_grid,lon_grid=np.meshgrid(lat,lon)\n",
    "on_land=globe.is_land(lat_grid, lon_grid).T\n",
    "\n",
    "coords=xr.where(on_land,data,np.nan) \n",
    "data=coords\n",
    "recon['spei_mean'].data=data.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2adc8ec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load neukom data, compute global mean surface temperature and core monsoon region temperature anomaly\n",
    "import glob\n",
    "paths=glob.glob('/home/mchoblet/neukom/*')\n",
    "lijst=[]\n",
    "for f in paths:\n",
    "    da=xr.open_dataset(f)\n",
    "    var=list(da.data_vars)[0]\n",
    "    da=da[var]\n",
    "    da=da.mean(axis=1)\n",
    "    #replace the time axis for having the same type as my reconstructiony\n",
    "    da['time']=recon['time'].values\n",
    "    lijst.append(da)\n",
    "    \n",
    "combined_dataset = xr.concat(lijst, dim='dataset')\n",
    "combined_dataset=(combined_dataset-combined_dataset.sel(time=slice('0850','1850')).mean('time')).mean('dataset')\n",
    "#combined_dataset=combined_dataset.mean('dataset')\n",
    "\n",
    "#globalmean surface temperature\n",
    "neu_glob=globalmean(combined_dataset)\n",
    "#regional mean\n",
    "neu_reg=regionalmean(combined_dataset)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d7f605d",
   "metadata": {},
   "source": [
    "# Video 1, climate fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8faf8034",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmaps={'tsurf_mean':'RdBu_r','prec_mean':'BrBG','spei_mean':'PuOr','d18O_mean':'PiYG_r'}\n",
    "units={'tsurf_mean':'$\\Delta$ [K]','prec_mean':'$\\Delta$ [mm/month]','spei_mean':'$\\Delta$ [1]','d18O_mean':'$\\Delta$ [‰]','pdsi_mean':'$\\Delta$ [1]'}\n",
    "var={'prec_mean':'Precipitation','spei_mean':'SPEI','tsurf_mean':'Temperature', 'd18O_mean':'$\\delta^{18}$O$_{prec}$'}\n",
    "maxs={'tsurf_mean':0.5,'prec_mean':5,'spei_mean':2.5,'d18O_mean':0.5}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "390772a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "resampled=recon.resample(time='25YS',label='right').mean()\n",
    "test=xr.cftime_range(start='0025',end='2000',freq='25YS')\n",
    "\n",
    "#resampled=recon.resample(time='50YS',label='right').mean()\n",
    "#test=xr.cftime_range(start='0050',end='2000',freq='50YS')\n",
    "\n",
    "resampled['time']=test.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "213561cc",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "map_extent=[250,340,-65,32]\n",
    "labelsize=20\n",
    "titlesize=23\n",
    "\n",
    "! rm -rf /home/mchoblet/paper_plots/videos/video1/*\n",
    "for i,t in enumerate(resampled.time):\n",
    "    fig, axes = plt.subplots(2, 2, figsize=(12,10), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "    year=t.dt.year.values.tolist()\n",
    "    string=str(year)+' CE'\n",
    "    string=f\"{string:>10}\"\n",
    "    fig.suptitle(string,fontsize=30,ha='center',y=0.96)\n",
    "    for iii,(ii,v) in enumerate(maxs.items()):\n",
    "        #ii,i,v    \n",
    "        ax=axes.flatten()[iii]\n",
    "        ax.set_extent(map_extent, crs=ccrs.PlateCarree())\n",
    "        ax.coastlines()\n",
    "        plot=resampled[ii].isel(time=i).plot(vmax=v,vmin=-v,ax=ax,cmap=cmaps[ii],add_colorbar=False,levels=21)\n",
    "        year=resampled.time.dt.year.values[i]\n",
    "        ax.set_ylabel('',fontsize=30,rotation=0)\n",
    "        ax.set_xlabel('',fontsize=30,rotation=0)\n",
    "        \n",
    "        p0=ax.get_position().get_points().flatten()\n",
    "        ax_cbar1 = fig.add_axes([p0[2]+0.03, p0[1], 0.01, p0[3]-p0[1]])\n",
    "        cb=fig.colorbar(plot, cax=ax_cbar1, orientation='vertical',cmap=cmaps[ii],extend='neither')\n",
    "        cb.ax.tick_params(labelsize=labelsize-3)\n",
    "\n",
    "        if 'temp' in ii.lower() or 'tsurf' in ii.lower() or 'd18' in ii.lower():\n",
    "            labelpad=-10\n",
    "        else: labelpad=2\n",
    "        \n",
    "        cb.set_label(units[ii],fontsize=labelsize,labelpad=labelpad)\n",
    "        cb.set_ticks([-v*4/5,0,v*4/5])\n",
    "        ax.set_title(ii,fontsize=20)\n",
    "        ax.set_title(var[ii],fontsize=titlesize)\n",
    "        fig.savefig(\"/home/mchoblet/paper_plots/revision/videos/video1/file%d.png\" % i, dpi=200,bbox_inches='tight')\n",
    "        plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05f6d3ab",
   "metadata": {},
   "source": [
    "## How to make video\n",
    "\n",
    "ffmpeg -framerate 2 -i file%02d.png -c:v libx264 -pix_fmt yuv420p video1.mp4\n",
    "\n",
    "Choose a slow framerate as reconstruction already resampled to 50 year means (40 images -> 20 seconds in total)\n",
    "\n",
    "ffmpeg -r 4 -i 'file%d.png' -vf \"scale=trunc(iw/2)*2:trunc(ih/2)*2\" -c:v libx264 -pix_fmt yuv420p  video1.mp4 \n",
    "(scaling needed for some x264 issue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6d7349d",
   "metadata": {},
   "source": [
    "# Video2: Monsoon Index for Precipitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb01bc66",
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_mean_filt=regionalmean(recon['prec_mean'])\n",
    "reg_mean_filt.data=quick_filt(reg_mean_filt,cutoff=1/50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4696e93c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#add regional temperature\n",
    "reg_mean_temp_filt=regionalmean(recon['tsurf_mean'])\n",
    "reg_mean_temp_filt.data=quick_filt(reg_mean_temp_filt,cutoff=1/50)\n",
    "\n",
    "#neukom\n",
    "neu_reg_filt=neu_reg.copy(deep=True)\n",
    "neu_reg_filt.data=quick_filt(neu_reg,cutoff=1/50)\n",
    "neu_glob_filt=neu_reg.copy(deep=True)\n",
    "neu_glob_filt.data=quick_filt(neu_glob,cutoff=1/50)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d65d06d7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from matplotlib.gridspec import GridSpec\n",
    "\n",
    "labels=[1,400,800,1200,1600,2000]\n",
    "map_extent=[250,340,-65,32]\n",
    "figsize=(20,10)\n",
    "lw=3\n",
    "\n",
    "data=resampled['prec_mean']\n",
    "\n",
    "for i,t in enumerate(data.time):\n",
    "    fig = plt.figure(figsize=figsize,constrained_layout=0)\n",
    "    gs = GridSpec(8,2, figure=fig)\n",
    "    ax = fig.add_subplot(gs[:,0],projection=ccrs.PlateCarree())\n",
    "\n",
    "    ax.add_feature(cfeature.COASTLINE,zorder=1)\n",
    "    gl=ax.axes.gridlines(color='black',alpha=0.5,linestyle='--',draw_labels=True, dms=True, x_inline=False, y_inline=False)\n",
    "    gl.top_labels = False\n",
    "    gl.right_labels = False\n",
    "\n",
    "    ax.axes.xaxis.set_visible(False)\n",
    "    ax.axes.yaxis.set_visible(False)\n",
    "    ax.axes.xaxis.set_ticks([])\n",
    "    ax.axes.yaxis.set_ticks([])\n",
    "    ax.set_extent(map_extent, crs=ccrs.PlateCarree())\n",
    "    ax.plot([285.5,312.5,312.5,285.5,285.5],[-17.5,-17.5,-5,-5,-17.5], transform=ccrs.PlateCarree(),color='blue',zorder=0)\n",
    "    ax.add_feature(cfeature.COASTLINE)\n",
    "    ax.axes.gridlines(color='black',alpha=0.5,linestyle='--')\n",
    "    gl.xlabel_style = {'size': 20}\n",
    "    gl.ylabel_style = {'size': 20}\n",
    "\n",
    "        \n",
    "    title='Precipitation'\n",
    "    vmax=6.2; vmin=-3\n",
    "    cmap='BrBG' #BrBG_r for precipitation\n",
    "    label=r'$\\Delta$Prec [mm/month]'\n",
    "    \n",
    "    sel=resampled['prec_mean'].isel(time=i)\n",
    "    plot=sel.plot(vmax=5,vmin=-5,ax=ax,cmap=cmap,add_colorbar=False)\n",
    "    year=sel.time.dt.year.values.tolist()\n",
    "\n",
    "    string=str(year)+' CE'\n",
    "    string=f\"{string:>10}\"\n",
    "    ax.set_title(string,fontsize=30)\n",
    "    #prox_i=proxies[i].dropna('site')\n",
    "    #scatter=ax.scatter(prox_i['lon'],prox_i['lat'],c=prox_i,s=size, cmap='PiYG_r', transform=ccrs.PlateCarree(), edgecolors='black',vmin=-2,vmax=2)\n",
    "\n",
    "    ax1 = fig.add_subplot(gs[1:5,1])        \n",
    "\n",
    "    #ax1.plot(time,data,label='LMR',lw=lw)\n",
    "    #if i==0:\n",
    "    lns1=ax1.plot(reg_mean_filt.values[:((i+1)*25)],color='black',lw=lw,label='$\\Delta$ Prec (recon)')\n",
    "    ax1.plot(reg_mean_filt.values,alpha=0)\n",
    "    ax1.set_xlim(-100,1900)\n",
    "    \n",
    "    \n",
    "    ax2 = ax1.twinx() \n",
    "    #usual linestyle is the best (else it looks tordue)\n",
    "    lns2=ax2.plot(reg_mean_temp_filt.values[:((i+1)*25)],lw=lw,ls='-',label='$\\Delta$ Temp (recon)')\n",
    "    lns3=ax2.plot(neu_reg_filt.values[:((i+1)*25)],lw=lw,ls='-',label='$\\Delta$ Temp (Neukom 2018)')\n",
    "    lns4=ax2.plot(neu_glob_filt.values[:((i+1)*25)],lw=lw,ls='-',label='$\\Delta$ GMST (Neukom 2018)')\n",
    "    lns = lns1+lns2+lns3+lns4\n",
    "    labs = [l.get_label() for l in lns]\n",
    "    \n",
    "    ax1.legend(lns, labs,bbox_to_anchor=(1.1, -0.27),ncols=2,fontsize=20,frameon=True,framealpha=1,edgecolor='black')\n",
    "    \n",
    "    xticks=np.arange(0,1,0.2)*2500\n",
    "    labels=[1,500,1000,1500,2000]\n",
    "    ax1.set_xticks(xticks,labels)\n",
    "    \n",
    "    p0=ax.get_position().get_points().flatten()\n",
    "    ax.plot([285.5,312.5,312.5,285.5,285.5],[-17.5,-17.5,-5,-5,-17.5], transform=ccrs.PlateCarree(),color='black',ls='--',lw=3)\n",
    "    \n",
    "    #ax1.set_xticklabels(labels,fontsize=25)\n",
    "    ax1.set_xlabel('Year (CE)',fontsize=30)\n",
    "    ax1.set_ylabel('Prec [mm/month]',fontsize=25)\n",
    "    ax2.set_ylabel('Temp [K]',fontsize=25)\n",
    "    \n",
    "    #hacky way to get the yaxis aligned somehow\n",
    "    ax1.set_yticks([-2,0,2,4,6])\n",
    "    ax2.set_yticks([-0.25,0,0.25,0.5,0.75])\n",
    "    #add a relative padding\n",
    "    ax1.set_ylim(-2-8/10,6+8/10)\n",
    "    ax2.set_ylim(-0.25-1/10,0.75+1/10)\n",
    "    \n",
    "    \n",
    "    ax2.yaxis.set_tick_params(labelsize=20)\n",
    "    ax1.yaxis.set_tick_params(labelsize=20)\n",
    "    ax2.spines['right'].set_visible(True)\n",
    "    #alignYaxes([ax1,ax2])\n",
    "    ax2.set_yticks(np.linspace(ax2.get_yticks()[0], ax2.get_yticks()[-1], len(ax1.get_yticks())))\n",
    "    ax2.grid(None)\n",
    "    \n",
    "    plt.title('Core monsoon region',fontsize=30,y=1.05)    \n",
    "    \n",
    "    ax_cbar1 = fig.add_axes([p0[0], p0[1]-0.08, (p0[2]-p0[0]), 0.02])\n",
    "    cb=fig.colorbar(plot, cax=ax_cbar1, orientation='horizontal',cmap=cmap)\n",
    "    cb.ax.tick_params(labelsize=20)\n",
    "    cb.set_ticks([-vmin*4/5,0,vmax*4/5])\n",
    "    label=r'$\\Delta$ Prec [mm/month]'\n",
    "    cb.set_label(label,fontsize=30)\n",
    "    vmax=5;vmin=-vmax\n",
    "    cb.set_ticks([vmin*4/5,0,vmax*4/5])\n",
    "    ax1.set_xticklabels(labels,fontsize=25)\n",
    "    \n",
    "    plt.savefig(\"/home/mchoblet/paper_plots/revision/videos/video2/file%d.png\" % i, dpi=200,bbox_inches='tight')\n",
    "    plt.show()\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1713315d",
   "metadata": {},
   "source": [
    "# Video3: Add Speleothem Anomalies\n",
    "\n",
    "* I need to load the speleothem data externally for getting the location of the proxy records, which were not stored during the reconstruction. But I use the resampled proxy record data (which is part of the reconstruction file).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a4e5c5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "speleos=xr.open_dataset(folder+files[0])['proxies_res_10']\n",
    "speleos=speleos.rename({'time_res_10':'time'})\n",
    "speleos=speleos-speleos.sel(time=slice('0850','1850')).mean('time')\n",
    "speleos=speleos.resample(time='25YS').mean()\n",
    "speleos['time']=test.values\n",
    "\n",
    "#have to add lats and lons as I forgot to add that in the PaleoDA code\n",
    "speleo_d=xr.open_dataset('/home/mchoblet/paleoda/data/proxy_dbs/paper/speleos.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47436eb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "array1=speleos['site_10'].values.tolist()\n",
    "array2=speleo_d['geo_siteName'].values.tolist()\n",
    "index_array = list(range(len(array2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a916a933",
   "metadata": {},
   "outputs": [],
   "source": [
    "index_array.sort(key=lambda x: array1.index(array2[x]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfe473b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "lats=speleo_d['geo_meanLat'][index_array]\n",
    "lons=speleo_d['geo_meanLon'][index_array]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "963084bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "speleo_d['geo_siteName'][index_array].values==speleos['site_10'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "320a7890",
   "metadata": {},
   "outputs": [],
   "source": [
    "speleos['lat']=('site_10',lats.values)\n",
    "speleos['lon']=('site_10',lons.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd47d568",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from matplotlib.gridspec import GridSpec\n",
    "\n",
    "labels=[1,400,800,1200,1600,2000]\n",
    "map_extent=[250,340,-65,32]\n",
    "figsize=(20,10)\n",
    "lw=3\n",
    "size=200\n",
    "\n",
    "data=resampled['prec_mean']\n",
    "\n",
    "for i,t in enumerate(data.time):\n",
    "    fig = plt.figure(figsize=figsize,constrained_layout=0)\n",
    "    gs = GridSpec(8,2, figure=fig)\n",
    "    ax = fig.add_subplot(gs[:,0],projection=ccrs.PlateCarree())\n",
    "\n",
    "    ax.add_feature(cfeature.COASTLINE,zorder=1)\n",
    "    gl=ax.axes.gridlines(color='black',alpha=0.5,linestyle='--',draw_labels=True, dms=True, x_inline=False, y_inline=False)\n",
    "    gl.top_labels = False\n",
    "    gl.right_labels = False\n",
    "\n",
    "    ax.axes.xaxis.set_visible(False)\n",
    "    ax.axes.yaxis.set_visible(False)\n",
    "    ax.axes.xaxis.set_ticks([])\n",
    "    ax.axes.yaxis.set_ticks([])\n",
    "    ax.set_extent(map_extent, crs=ccrs.PlateCarree())\n",
    "    ax.plot([285.5,312.5,312.5,285.5,285.5],[-17.5,-17.5,-5,-5,-17.5], transform=ccrs.PlateCarree(),color='blue',zorder=0)\n",
    "    ax.add_feature(cfeature.COASTLINE)\n",
    "    ax.axes.gridlines(color='black',alpha=0.5,linestyle='--')\n",
    "    gl.xlabel_style = {'size': 20}\n",
    "    gl.ylabel_style = {'size': 20}\n",
    "\n",
    "        \n",
    "    title='Precipitation'\n",
    "    vmax=6.2; vmin=-3\n",
    "    cmap='BrBG' #BrBG_r for precipitation\n",
    "    label=r'$\\Delta$Prec [mm/month]'\n",
    "    \n",
    "    sel=resampled['prec_mean'].isel(time=i)\n",
    "    plot=sel.plot(vmax=5,vmin=-5,ax=ax,cmap=cmap,add_colorbar=False)\n",
    "    year=sel.time.dt.year.values.tolist()\n",
    "\n",
    "    string=str(year)+' CE'\n",
    "    string=f\"{string:>10}\"\n",
    "    ax.set_title(string,fontsize=30)\n",
    "    \n",
    "    prox_i=speleos.isel(time=i).dropna('site_10')\n",
    "    vmax_18=1\n",
    "    scatter=ax.scatter(prox_i['lon'],prox_i['lat'],c=prox_i,s=size, cmap='PiYG_r', transform=ccrs.PlateCarree(), edgecolors='black',vmin=-vmax_18,vmax=vmax_18)\n",
    "    \n",
    "    ax1 = fig.add_subplot(gs[1:5,1])        \n",
    "\n",
    "    #ax1.plot(time,data,label='LMR',lw=lw)\n",
    "    #if i==0:\n",
    "    lns1=ax1.plot(reg_mean_filt.values[:((i+1)*25)],color='black',lw=lw,label='$\\Delta$ Prec (recon)')\n",
    "    ax1.plot(reg_mean_filt.values,alpha=0)\n",
    "    ax1.set_xlim(-100,1900)\n",
    "    \n",
    "    \n",
    "    ax2 = ax1.twinx() \n",
    "    #usual linestyle is the best (else it looks tordue)\n",
    "    lns2=ax2.plot(reg_mean_temp_filt.values[:((i+1)*25)],lw=lw,ls='-',label='$\\Delta$ Temp (recon)')\n",
    "    lns3=ax2.plot(neu_reg_filt.values[:((i+1)*25)],lw=lw,ls='-',label='$\\Delta$ Temp (Neukom 2018)')\n",
    "    lns4=ax2.plot(neu_glob_filt.values[:((i+1)*25)],lw=lw,ls='-',label='$\\Delta$ GMST (Neukom 2018)')\n",
    "    lns = lns1+lns2+lns3+lns4\n",
    "    labs = [l.get_label() for l in lns]\n",
    "    \n",
    "    ax1.legend(lns, labs,bbox_to_anchor=(1.1, -0.27),ncols=2,fontsize=20,frameon=True,framealpha=1,edgecolor='black')\n",
    "    \n",
    "    xticks=np.arange(0,1,0.2)*2500\n",
    "    labels=[1,500,1000,1500,2000]\n",
    "    ax1.set_xticks(xticks,labels)\n",
    "    \n",
    "    p0=ax.get_position().get_points().flatten()\n",
    "    ax.plot([285.5,312.5,312.5,285.5,285.5],[-17.5,-17.5,-5,-5,-17.5], transform=ccrs.PlateCarree(),color='black',ls='--',lw=3)\n",
    "    \n",
    "    #ax1.set_xticklabels(labels,fontsize=25)\n",
    "    ax1.set_xlabel('Year (CE)',fontsize=30)\n",
    "    ax1.set_ylabel('Prec [mm/month]',fontsize=25)\n",
    "    ax2.set_ylabel('Temp [K]',fontsize=25)\n",
    "    \n",
    "    #hacky way to get the yaxis aligned somehow\n",
    "    ax1.set_yticks([-2,0,2,4,6])\n",
    "    ax2.set_yticks([-0.25,0,0.25,0.5,0.75])\n",
    "    #add a relative padding\n",
    "    ax1.set_ylim(-2-8/10,6+8/10)\n",
    "    ax2.set_ylim(-0.25-1/10,0.75+1/10)\n",
    "    \n",
    "    \n",
    "    ax2.yaxis.set_tick_params(labelsize=20)\n",
    "    ax1.yaxis.set_tick_params(labelsize=20)\n",
    "    ax2.spines['right'].set_visible(True)\n",
    "    #alignYaxes([ax1,ax2])\n",
    "    ax2.set_yticks(np.linspace(ax2.get_yticks()[0], ax2.get_yticks()[-1], len(ax1.get_yticks())))\n",
    "    ax2.grid(None)\n",
    "    \n",
    "    plt.title('Core monsoon region',fontsize=30,y=1.05)    \n",
    "    \n",
    "    ax_cbar1 = fig.add_axes([p0[0], p0[1]-0.08, (p0[2]-p0[0]), 0.02])\n",
    "    cb=fig.colorbar(plot, cax=ax_cbar1, orientation='horizontal',cmap=cmap)\n",
    "    cb.ax.tick_params(labelsize=20)\n",
    "    label=r'$\\Delta$ Prec [mm/month]'\n",
    "    cb.set_label(label,fontsize=30)\n",
    "    vmax=5;vmin=-vmax\n",
    "    cb.set_ticks([vmin*4/5,0,vmax*4/5])\n",
    "    ax1.set_xticklabels(labels,fontsize=25)\n",
    "    \n",
    "    ax_cbar1 = fig.add_axes([p0[0], p0[1]-0.2, (p0[2]-p0[0]), 0.02])\n",
    "    cb=fig.colorbar(scatter, cax=ax_cbar1, orientation='horizontal')\n",
    "    cb.ax.tick_params(labelsize=20)\n",
    "    cb.set_label(u'Speleothem $\\Delta(\\delta^{18}$O) [‰]',fontsize=30)\n",
    "    cb.ax.invert_xaxis()\n",
    "    vmax=vmax_18;vmin=-vmax\n",
    "    cb.set_ticks([vmin*4/5,0,vmax*4/5])\n",
    "    ax1.set_xticklabels(labels,fontsize=25)\n",
    "    \n",
    "    plt.savefig(\"/home/mchoblet/paper_plots/revision/videos/video3/file%d.png\" % i, dpi=200,bbox_inches='tight')\n",
    "    plt.show()\n",
    "    plt.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "paleoda_kernel",
   "language": "python",
   "name": "paleoda_kernel"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}