notebooks/calculate_global_statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4b762609-8a73-40a3-82d1-638bdc2a4e36",
   "metadata": {},
   "source": [
    "# Correlation statistics for global variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fee43bc-ddc5-4243-8e93-41823353767c",
   "metadata": {},
   "source": [
    "This notebook calculates slope, p value, standard deviation of the slope, and correlation coefficients for the linear and exponential fit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf747a04-a34c-432c-bf97-765a85db0232",
   "metadata": {
    "tags": []
   },
   "source": [
    "Loading libraires:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2aed15-88fd-4cb8-ac97-1518c3b28f68",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import netCDF4 as nc\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy import stats\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c013d05-3876-4465-8e26-9b6586667e9b",
   "metadata": {},
   "source": [
    "Define functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90de2d0d-81a9-4640-8fc9-71967662b1af",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def read_txt(fname, ncolumns, dtype):\n",
    "    fdata   = open(fname)\n",
    "    flines  = np.array(fdata.readlines())\n",
    "    data_in = np.empty([len(flines)])\n",
    "    if dtype == 'string':\n",
    "        data_in = np.empty([len(flines)], dtype=object)\n",
    "    for i in np.arange(len(flines)):\n",
    "        data_in[i]=flines[i].split()[ncolumns-1]      \n",
    "    return(data_in)\n",
    "    fdata.close()\n",
    "    \n",
    "def read_data(fname, ncolumns):\n",
    "    if ncolumns > 1:\n",
    "        fdata   = open(fname)\n",
    "        flines  = np.array(fdata.readlines())\n",
    "        data_in = np.empty([len(flines), ncolumns])\n",
    "        for i in np.arange(len(flines)):\n",
    "            data_in[i,:]=flines[i].split()\n",
    "    else:\n",
    "        with open(fname) as flist:\n",
    "            data_in = np.array(flist.readlines())        \n",
    "    return(data_in)\n",
    "    fdata.close()\n",
    "\n",
    "def read_nc(fname, varname, ifmerged):\n",
    "    if ifmerged:\n",
    "        f  = nc.MFDataset(fname)\n",
    "    else:\n",
    "        f  = nc.Dataset(fname)\n",
    "    var    = f.variables[varname]\n",
    "    lon    = f.variables['lon']\n",
    "    lat    = f.variables['lat']\n",
    "    time   = f.variables['time']\n",
    "    dates  = nc.num2date(time[:], time.units, time.calendar)\n",
    "    return var, lon, lat, time, dates\n",
    "\n",
    "def read_mask(fname, varname):\n",
    "    f = nc.Dataset(fname)\n",
    "    var  = f.variables[varname][:]\n",
    "    return var"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ded0b8a-fd13-4a73-9014-f526e3577026",
   "metadata": {},
   "source": [
    "Define paths and enviromental variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d52f0475-3a6d-44d6-be60-fcfc942260f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "resolution   = 'low_res' \n",
    "seasons      = ['ANNUAL', 'DJF', 'MAM', 'JJA', 'SON'] #'SON' #'SON' #None #\n",
    "member_names = f'models_{resolution}.txt'\n",
    "varname      = 'tos'\n",
    "ssps         = ['ssp126', 'ssp245', 'ssp370', 'ssp585']\n",
    "avg_period   = 10 #, 20, 30]\n",
    "ref_ntimes   = 51\n",
    "nssps = len(ssps)\n",
    "ystart_to_plot  = 2010\n",
    "ystart = int(ystart_to_plot - 1950)\n",
    "\n",
    "root = '..'\n",
    "path2data = f'{root}/data/data_txt/'  \n",
    "path2info  = f'{root}/data/data_info/' "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd92982-5ef7-4ea4-93d6-3b8cc3f294bc",
   "metadata": {
    "tags": []
   },
   "source": [
    "Read list of the models to be processed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa23e116-0528-48c2-85df-0295dbd0c7cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "member_list = f'{path2info}/{member_names}'\n",
    "members  = read_txt(member_list, 1, 'string')\n",
    "nmembers = len(members)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b850ca80-86e7-4da6-8460-f997178ec136",
   "metadata": {},
   "source": [
    "Main code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "792af0a8-47a6-48c2-ac7b-0a32c2a7cd06",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Loop over seasons\n",
    "for season in seasons:       \n",
    "    print(f'Working on the season {season}')\n",
    "    \n",
    "    Beta   = np.empty((nssps),dtype=float)\n",
    "    R2_lin = np.empty((nssps),dtype=float)\n",
    "    R2_exp = np.empty((nssps),dtype=float)\n",
    "    p_value = np.empty((nssps),dtype=float)\n",
    "    std_value = np.empty((nssps),dtype=float)       \n",
    "    \n",
    "        \n",
    "    for s in np.arange(nssps):\n",
    "        ssp = ssps[s]\n",
    "        print(\"        Working on the ssp: \" + ssp)      \n",
    "\n",
    "        # Read reference data = global mean temperature\n",
    "        ref_fname = f'{path2data}/{resolution}/global/tas/tas_{ssp}_ANNUAL.csv'\n",
    "        ref_data = read_data(ref_fname, nmembers)\n",
    "\n",
    "        # Read model data\n",
    "        dirpath  = f'{resolution}/global/{varname}/' \n",
    "        filename = f'{varname}_{ssp}_{season}'\n",
    "        infile   = f'{path2data}/{dirpath}/{filename}.csv'\n",
    "        data     = read_data(infile, nmembers)         \n",
    "\n",
    "        # Linear fit\n",
    "        x = ref_data[ystart::,:].reshape(-1, 1)\n",
    "        y = data[ystart::,:].reshape(-1, 1)\n",
    "        Beta[s], intercept, rvalue, p_value[s], std_value[s] = stats.linregress(x[:,0], y[:,0])\n",
    "        R2_lin[s] = rvalue*rvalue  \n",
    "\n",
    "        # Exponential fit\n",
    "        def func(x, a, b): return a * np.exp(b * x)\n",
    "        popt, pcov = curve_fit(func, x[:,0], y[:,0])\n",
    "        x_fitted = np.linspace(np.min(x), np.max(x), 150)\n",
    "        y_fitted = func(x_fitted, *popt)         \n",
    "        residuals = y[:,0] - func(x[:,0], *popt)\n",
    "        res = np.sum(residuals**2)\n",
    "        tot = np.sum((y[:,0]-np.mean(y[:,0]))**2)\n",
    "        R2_exp[s] = 1 - (res / tot)\n",
    "\n",
    "    # Save statistical info in a txt file\n",
    "    print(\"    Saving the data in the output file.\")    \n",
    "    path2out  = f'{path2data}/{resolution}/global/{varname}'\n",
    "    if not os.path.exists(path2out):\n",
    "        os.makedirs(path2out)\n",
    "        \n",
    "    print(Beta)\n",
    "\n",
    "    np.savetxt(f'{path2out}/Beta_{varname}_{season}.csv',      Beta,      delimiter =\" \", fmt='%6.2f')\n",
    "    np.savetxt(f'{path2out}/R2lin_{varname}_{season}.csv',     R2_lin,    delimiter =\" \", fmt='%6.2f')\n",
    "    np.savetxt(f'{path2out}/R2exp_{varname}_{season}.csv',     R2_exp,    delimiter =\" \", fmt='%6.2f')\n",
    "    np.savetxt(f'{path2out}/p_value_{varname}_{season}.csv',   p_value,   delimiter =\" \", fmt='%6.5f')\n",
    "    np.savetxt(f'{path2out}/std_value_{varname}_{season}.csv', std_value, delimiter =\" \", fmt='%6.5f')        "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}