{ "cells": [ { "cell_type": "markdown", "id": "60dc9908-b54e-4702-a2ac-7b2a39ceb78f", "metadata": {}, "source": [ "## *NOTE* Frequency indices T?90, R??mm, N??? have errors in them. Output files deleted." ] }, { "cell_type": "markdown", "id": "4176fb9d-fa11-4c3f-9fe3-ab8f85c8a2c6", "metadata": {}, "source": [ "Calculate the climatology and anomalies for each index. \n", "Then calculate the statistical significance to plot results\n", "\n", "Using ARISE-SAI runs with 5% seeding. Because only 5 runs of SSP245 have tmax and tmin, subtract the ensemble mean climatology from all temperatures." ] }, { "cell_type": "code", "execution_count": 2, "id": "fc0a7e24-8f74-481e-ad77-a5a92941f37f", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "code", "execution_count": 3, "id": "8716581d-05e0-4263-938e-508dbf434365", "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy as np # data arrays\n", "import xarray as xr # data array manipulation\n", "import pandas as pd\n", "import os\n", "import glob # to enable pattern matching\n", "from scipy import stats" ] }, { "cell_type": "code", "execution_count": 4, "id": "f911edaf-07da-4e30-89e0-55b9f6849768", "metadata": {}, "outputs": [], "source": [ "iDir = \"/glade/work/maritye/Data/ARISE-MCB/ETCCDI/\"\n", "inSSP = '/glade/work/maritye/Data/ARISE-SAI/ETCCDI/SSP245/'\n", "oDir = '/glade/work/maritye/Data/ARISE-MCB/ETCCDI/Anomalies/'\n", "\n", "varnam = ['PRECT', 'TREFHTMX', 'TREFHTMN']" ] }, { "cell_type": "markdown", "id": "934c4a5e-2606-46ed-9544-cb8630a14365", "metadata": {}, "source": [ "Can't quite work out how to loop the directories and the separate lists of indices. So just do three repeated calculations for anomalies and significances" ] }, { "cell_type": "markdown", "id": "c58ec6c6-44de-486d-8027-081e24e3bee3", "metadata": {}, "source": [ "Create standard information to include with each data set" ] }, { "cell_type": "code", "execution_count": 77, "id": "a6652955-d0ba-4821-928e-dc2149e181c7", "metadata": {}, "outputs": [], "source": [ "years = np.arange(2035,2070)\n", "dates= '2035-2069'\n", "member=range(10)\n", "lat = np.linspace(-90,90,num=192)\n", "lon = np.linspace(0,358.75,num=288)\n", "dims = ('member', 'lat', 'lon')\n", "coords = dict(member = member, lat=lat, lon=lon)\n", "attribs = dict(description='Extreme Index Anomalies based on ETCCDI definitions. MCB scenarios with 5% global seeding', \n", " history='Created by Mari Tye November 2023.' ),\n" ] }, { "cell_type": "code", "execution_count": 108, "id": "67f7f85b-812d-4089-bb6a-1591706a341b", "metadata": {}, "outputs": [], "source": [ "pindices = ['CDD','CWD','RX1D','N95','N99','R20mm','RX5D','PRCPTOT', 'SDII','R10mm'] # something funny with P95 and P99 ignore these 2 for MCB runs\n", "txindices = ['TXX', 'TXN', 'TX10', 'TX90'] #calculated the fixed threshold indices but they aren't very meaningful at a global scale so not including in paper or significance tests.\n", "tnindices = ['TNX', 'TNN', 'TN10', 'TN90']" ] }, { "cell_type": "markdown", "id": "9668dfbd-f2d4-458b-a738-990ea8e6640f", "metadata": {}, "source": [ "Precipitation anomalies and significance of changes" ] }, { "cell_type": "code", "execution_count": 65, "id": "9f339d23-5865-48fd-8cff-185932a86add", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CDD\n", "CWD\n", "RX1D\n", "N95\n", "N99\n", "R20mm\n", "RX5D\n", "PRCPTOT\n", "SDII\n", "R10mm\n" ] } ], "source": [ "for i in range(len(pindices)):\n", " print(pindices[i])\n", " \n", " # this creates a vector, but only want first element to use as file name\n", " praw_ssp = glob.glob(inSSP + varnam[0] + '/*' + pindices[i] + '*.nc')[0]\n", " praw_mcb = glob.glob(iDir + varnam[0] + '/*' + pindices[i] + '*.nc')[0]\n", " \n", " with xr.open_dataset(praw_ssp) as ds_ssp:\n", " cont = ds_ssp.sel(year=slice('2019-01-01','2039-12-31'))\n", " ssp_ec = ds_ssp.sel(year=slice('2049-01-01', '2069-12-31'))\n", " clim = cont.mean(dim='year')\n", " \n", " with xr.open_dataset(praw_mcb) as ds_mcb:\n", " mcb_ec = ds_mcb.sel(year=slice('2049-01-01', '2069-12-31'))\n", " \n", " \n", " \n", " bc_anom = cont - clim\n", " ssp_anom = ssp_ec - clim\n", " mcb_anom = mcb_ec - clim\n", " \n", " \n", " sigs = xr.DataArray(None, coords=dict(member=member, lat = lat, lon = lon), dims=dims, name='pval_ssp')\n", " sigm = xr.DataArray(None, coords=dict(member=member, lat = lat, lon = lon), dims=dims, name='pval_mcb')\n", "\n", " for em in range(10):\n", " \n", " index_c = bc_anom[pindices[i]].sel(member = em)\n", " index_e = ssp_anom[pindices[i]].sel(member = em)\n", " index_m = mcb_anom[pindices[i]].sel(member = em)\n", " foo = stats.ttest_ind(index_c, index_e).pvalue\n", " boo = stats.ttest_ind(index_c, index_m).pvalue\n", " sigs[em,:,:,] = xr.DataArray(foo)\n", " sigm[em,:,:,] = xr.DataArray(boo)\n", " \n", " # create a merged file of end of century ensemble mean anomalies and significances\n", " ssp_ec_mean = (ssp_anom.mean(dim=('year', 'member'))).to_array(name = f\"SSP_{str(pindices[i])}\")\n", " mcb_ec_mean = (mcb_anom.mean(dim=('year', 'member'))).to_array(name = f\"MCB_{str(pindices[i])}\")\n", " ssp_sig = xr.DataArray(sigs.mean(dim='member'))\n", " mcb_sig = xr.DataArray(sigm.mean(dim='member'))\n", " \n", " anoms_sigs = xr.merge([ssp_ec_mean, mcb_ec_mean, ssp_sig, mcb_sig])\n", " anoms_sigs.assign_attrs(purpose = \"Ensemble mean climatology for producing figures\")\n", " flnm = iDir + 'Anomalies_Significances_2060_' + pindices[i] + '.nc'\n", " anoms_sigs.to_netcdf(flnm)\n", " \n", " \"\"\"\n", "# Save the anomaly time series for future reference.\n", " sspnm = praw_ssp.split(pindices[i])[0] + pindices[i] + '_2050-2069_anoms.nc'\n", " contnm = praw_ssp.split(pindices[i])[0] + pindices[i] + '_2020-2039_anoms.nc'\n", " mcbnm = praw_mcb.split(pindices[i])[0] + pindices[i] + '_2050-2069_anoms.nc'\n", " \n", " ssp_anom.to_netcdf(sspnm)\n", " mcb_anom.to_netcdf(mcbnm)\n", " bc_anom.to_netcdf(contnm)\n", "\"\"\" \n", "\n", " " ] }, { "cell_type": "code", "execution_count": 110, "id": "ebc3f5a8-0d4e-495d-aa0c-21aa466d54fa", "metadata": { "tags": [] }, "outputs": [], "source": [ "del(sigs)\n", "del(sigm)" ] }, { "cell_type": "code", "execution_count": 120, "id": "f0a98dd8-44a1-44af-818a-3ac52e3f3214", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TXX\n", "TXN\n", "TX10\n", "TX90\n" ] } ], "source": [ "for i in range(len(txindices)):\n", " print(txindices[i])\n", " \n", " # this creates a vector, but only want first element to use as file name\n", " praw_ssp = glob.glob(inSSP + varnam[1] + '/*' + txindices[i] + '*.nc')[0]\n", " praw_mcb = glob.glob(iDir + varnam[1] + '/*' + txindices[i] + '*.nc')[0]\n", " \n", " with xr.open_dataset(praw_ssp) as ds_ssp:\n", " ds_ssp = ds_ssp.rename({'members':'member'})\n", " cont = ds_ssp.sel(year=slice('2019-01-01','2039-12-31'))\n", " ssp_ec = ds_ssp.sel(year=slice('2049-01-01', '2069-12-31'))\n", " clim = cont.mean(dim='year')\n", " \n", " with xr.open_dataset(praw_mcb) as ds_mcb:\n", " mcb_ec = ds_mcb.sel(year=slice('2049-01-01', '2069-12-31'))\n", " \n", " \n", " \n", " bc_anom = cont - clim\n", " ssp_anom = ssp_ec - clim\n", " mcb_anom = mcb_ec - clim\n", " \n", " \n", " sigs = xr.DataArray(None, coords=dict(member=member, lat = lat, lon = lon), dims=dims, name='pval_ssp')\n", " sigm = xr.DataArray(None, coords=dict(member=member, lat = lat, lon = lon), dims=dims, name='pval_mcb')\n", "\n", " for em in range(10):\n", " \n", " index_c = bc_anom[txindices[i]].sel(member = em)\n", " index_e = ssp_anom[txindices[i]].sel(member = em)\n", " index_m = mcb_anom[txindices[i]].sel(member = em)\n", " foo = stats.ttest_ind(index_c, index_e).pvalue\n", " boo = stats.ttest_ind(index_c, index_m).pvalue\n", " sigs[em,:,:,] = xr.DataArray(foo)\n", " sigm[em,:,:,] = xr.DataArray(boo)\n", " \n", " # create a merged file of end of century ensemble mean anomalies and significances\n", " ssp_ec_mean = (ssp_anom.mean(dim=('year', 'member'))).to_array(name = f\"SSP_{str(txindices[i])}\")\n", " mcb_ec_mean = (mcb_anom.mean(dim=('year', 'member'))).to_array(name = f\"MCB_{str(txindices[i])}\")\n", " ssp_sig = xr.DataArray(sigs.mean(dim='member'))\n", " mcb_sig = xr.DataArray(sigm.mean(dim='member'))\n", " \n", " anoms_sigs = xr.merge([ssp_ec_mean, mcb_ec_mean, ssp_sig, mcb_sig])\n", " anoms_sigs.assign_attrs(purpose = \"Ensemble mean climatology for producing figures\")\n", " flnm = iDir + 'Anomalies_Significances_2060_' + txindices[i] + '.nc'\n", " anoms_sigs.to_netcdf(flnm)\n", " \n", "\n", "\n", " " ] }, { "cell_type": "markdown", "id": "8271e70d-6a18-46b2-9e38-bd2bc8cdf859", "metadata": { "tags": [] }, "source": [ "Daily minima Again using only the max, min and 90th/10th percentiles" ] }, { "cell_type": "code", "execution_count": 114, "id": "47a50f2f-7a2a-4a29-931b-f590080696b9", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TNX\n", "TNN\n", "TN10\n", "TN90\n" ] } ], "source": [ "for i in range(len(tnindices)):\n", " print(tnindices[i])\n", " \n", " # this creates a vector, but only want first element to use as file name\n", " praw_ssp = glob.glob(inSSP + varnam[2] + '/*' + tnindices[i] + '*.nc')[0]\n", " praw_mcb = glob.glob(iDir + varnam[2] + '/*' + tnindices[i] + '*.nc')[0]\n", " \n", " with xr.open_dataset(praw_ssp) as ds_ssp:\n", " ds_ssp = ds_ssp.rename({'members':'member'})\n", " cont = ds_ssp.sel(year=slice('2019-01-01','2039-12-31'))\n", " ssp_ec = ds_ssp.sel(year=slice('2049-01-01', '2069-12-31'))\n", " clim = cont.mean(dim='year')\n", " \n", " with xr.open_dataset(praw_mcb) as ds_mcb:\n", " mcb_ec = ds_mcb.sel(year=slice('2049-01-01', '2069-12-31'))\n", " \n", " \n", " \n", " bc_anom = cont - clim\n", " ssp_anom = ssp_ec - clim\n", " mcb_anom = mcb_ec - clim\n", " \n", " \n", " sigs = xr.DataArray(None, coords=dict(member=member, lat = lat, lon = lon), dims=dims, name='pval_ssp')\n", " sigm = xr.DataArray(None, coords=dict(member=member, lat = lat, lon = lon), dims=dims, name='pval_mcb')\n", "\n", " for em in range(10):\n", " \n", " index_c = bc_anom[tnindices[i]].sel(member = em)\n", " index_e = ssp_anom[tnindices[i]].sel(member = em)\n", " index_m = mcb_anom[tnindices[i]].sel(member = em)\n", " foo = stats.ttest_ind(index_c, index_e).pvalue\n", " boo = stats.ttest_ind(index_c, index_m).pvalue\n", " sigs[em,:,:,] = xr.DataArray(foo)\n", " sigm[em,:,:,] = xr.DataArray(boo)\n", " \n", " # create a merged file of end of century ensemble mean anomalies and significances\n", " ssp_ec_mean = (ssp_anom.mean(dim=('year', 'member'))).to_array(name = f\"SSP_{str(txindices[i])}\")\n", " mcb_ec_mean = (mcb_anom.mean(dim=('year', 'member'))).to_array(name = f\"MCB_{str(txindices[i])}\")\n", " ssp_sig = xr.DataArray(sigs.mean(dim='member'))\n", " mcb_sig = xr.DataArray(sigm.mean(dim='member'))\n", " \n", " anoms_sigs = xr.merge([ssp_ec_mean, mcb_ec_mean, ssp_sig, mcb_sig])\n", " anoms_sigs.assign_attrs(purpose = \"Ensemble mean climatology for producing figures\")\n", " flnm = iDir + 'Anomalies_Significances_2060_' + tnindices[i] + '.nc'\n", " anoms_sigs.to_netcdf(flnm)\n", " \n", "\n", "\n", " " ] }, { "cell_type": "code", "execution_count": 115, "id": "1b7167f2-6027-4e6b-a545-7d17c8c48945", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:   (lat: 192, lon: 288, variable: 1, members: 10)\n",
       "Coordinates:\n",
       "  * lat       (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0\n",
       "  * lon       (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8\n",
       "  * variable  (variable) object 'TN90'\n",
       "  * members   (members) int64 0 1 2 3 4 5 6 7 8 9\n",
       "Data variables:\n",
       "    SSP_TX90  (variable, lat, lon) float64 11.83 11.83 11.83 ... 5.09 5.075 5.12\n",
       "    MCB_TX90  (variable, lat, lon, members) float64 1.355e+03 ... 885.0\n",
       "    pval_ssp  (lat, lon) float64 0.01352 0.01352 0.01352 ... 0.1332 0.1292\n",
       "    pval_mcb  (lat, lon) float64 1.409e-76 1.409e-76 ... 4.484e-74 6.86e-74
" ], "text/plain": [ "\n", "Dimensions: (lat: 192, lon: 288, variable: 1, members: 10)\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 88.12 89.06 90.0\n", " * lon (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8\n", " * variable (variable) object 'TN90'\n", " * members (members) int64 0 1 2 3 4 5 6 7 8 9\n", "Data variables:\n", " SSP_TX90 (variable, lat, lon) float64 11.83 11.83 11.83 ... 5.09 5.075 5.12\n", " MCB_TX90 (variable, lat, lon, members) float64 1.355e+03 ... 885.0\n", " pval_ssp (lat, lon) float64 0.01352 0.01352 0.01352 ... 0.1332 0.1292\n", " pval_mcb (lat, lon) float64 1.409e-76 1.409e-76 ... 4.484e-74 6.86e-74" ] }, "execution_count": 115, "metadata": {}, "output_type": "execute_result" } ], "source": [ "anoms_sigs" ] }, { "cell_type": "code", "execution_count": null, "id": "8e6a0320-fa73-4ed0-9b62-8c62fcee4b24", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:miniconda-npl2023a_mrt]", "language": "python", "name": "conda-env-miniconda-npl2023a_mrt-py" }, "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.15" } }, "nbformat": 4, "nbformat_minor": 5 }