ERA5_HWMI_VT.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "orange-nerve",
   "metadata": {},
   "outputs": [],
   "source": [
    "#This script computes the HWMI based on daily fields for a specific period in cross-validation. \n",
    "import matplotlib\n",
    "from scipy import signal\n",
    "import os\n",
    "import datetime\n",
    "import time\n",
    "import copy\n",
    "import shutil\n",
    "import sys\n",
    "from function_read import *\n",
    "import numpy as np\n",
    "from scipy.stats import gaussian_kde\n",
    "import xarray as xr\n",
    "import calendar\n",
    "import locale\n",
    "from netCDF4 import num2date, date2num\n",
    "from joblib import Parallel, delayed\n",
    "import joblib\n",
    "from netCDF4 import num2date, date2num, Dataset\n",
    "from numpy import dtype\n",
    "from HWMI import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fe7963c",
   "metadata": {},
   "outputs": [],
   "source": [
    "var='tmax'\n",
    "expname='ERA5'\n",
    "memb_str='0' # there are not members for ERA5\n",
    "nrealisation=1\n",
    "# Region #\n",
    "reg_name ='Europe'\n",
    "min_lon = -15\n",
    "min_lat = 25\n",
    "max_lon = 60\n",
    "max_lat = 70\n",
    "lats_bnds = np.array([min_lat,max_lat])\n",
    "lons_bnds = np.array([min_lon, max_lon])\n",
    "# Season #\n",
    "season='15MJJA'\n",
    "season_start_day=[5,15]\n",
    "season_start_day=[8,31]\n",
    "nday=109 #15 of May to 31st of Aug\n",
    "# Period #\n",
    "ref_year1=1981\n",
    "ref_year2=1983\n",
    "nyear=len(range(ref_year1,ref_year2))+1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af59914d",
   "metadata": {},
   "outputs": [],
   "source": [
    "low=False # if data with 1 degree resolution is used (True) or not (False)\n",
    "\n",
    "if (low==True):\n",
    "    path=\"/work/csp/vt17420/ERA5/daily/russo/\"+var+\"_r360x180/\"\n",
    "    # Where the outputs will be saved\n",
    "    dirout=\"/work/csp/vt17420/ERA5/daily/russo/\"+var+\"_r360x180/HWMI/\"\n",
    "    lons_reg=np.arange(min_lon,max_lon,1)\n",
    "    lats_reg=np.arange(min_lat,max_lat,1)\n",
    "    nlon=len(lons_reg)\n",
    "    nlat=len(lats_reg)\n",
    "else:\n",
    "    path=\"/work/csp/vt17420/ERA5/daily/global/\"+var+\"/\"\n",
    "    dirout=\"/data/csp/vt17420/CLINT_PRODUCTS/\"+expname+\"_HWMI_\"+var+\"/\"\n",
    "    lons_reg=np.arange(min_lon,max_lon+0.25,0.25)\n",
    "    lats_reg=np.arange(min_lat,max_lat+0.25,0.25)\n",
    "    nlon=len(lons_reg)\n",
    "    nlat=len(lats_reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "formal-warrior",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=np.zeros((nlon,nlat,nyear,nday,1)) # daily data for all the years\n",
    "for iyear,year in enumerate(range(ref_year1,ref_year2)):\n",
    "    days_may=np.linspace(15, 31, num=17)\n",
    "    obs1=xr.open_dataset(path+var+\"_\"+str(year)+'05.nc')\n",
    "    if low==False:\n",
    "        obs1 = obs1.sel(latitude=slice(max_lat,min_lat),longitude=slice(min_lon,max_lon))\n",
    "    obs=obs1.sel(time=obs1.time.time.dt.day.isin(days_may))\n",
    "    # selection of the days in june, july and august\n",
    "    for mon in range(6,9):\n",
    "        obs1=xr.open_dataset(path+var+\"_\"+str(year)+str(mon).zfill(2)+'.nc')\n",
    "        if low==False:\n",
    "            obs1 = obs1.sel(latitude=slice(max_lat,min_lat),longitude=slice(min_lon,max_lon))\n",
    "        obs=xr.merge(([obs,obs1]))\n",
    "    data[:,:,iyear,:,0]=np.transpose(obs.to_array()[0,:,:,:],[2,1,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "innovative-swift",
   "metadata": {},
   "outputs": [],
   "source": [
    "### PERCENTILE THRESHOLD\n",
    "#percent_thresh = 95\n",
    "percent_thresh = 90\n",
    "cv='CV'\n",
    "cv_str = cv\n",
    "### MINIMAL DURATION OF A HW\n",
    "#duration_min = 5\n",
    "duration_min = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58c66a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "nmemb=data.shape[4]\n",
    "ndayseas = nday//duration_min +1\n",
    "HWMI = np.zeros((nyear,nmemb,nlat,nlon))\n",
    "ndayexedthreshold = np.zeros((nyear,nmemb,nlat,nlon))\n",
    "subHWarray = np.zeros((nyear, ndayseas, nmemb, nlat, nlon))\n",
    "#HW = np.zeros((nyear,nmemb,nday,nlat,nlon))\n",
    "#DDthreshold = np.zeros((nyear,nmemb,nlat,nlon))\n",
    "#fitsubHWarray = np.zeros((nyear, ndayseas, nmemb, nlat, nlon))\n",
    "#impossible_fit_list = []\n",
    "#sstMeanarray = np.zeros((nyear,nmemb,nlat,nlon))\n",
    "\n",
    "def parallelized_HWMIs_computation(ilat, ilon):    #, HWMI, ndayexedthreshold, DDthreshold, fitsubHWarray, subHWarray):\n",
    "        HWMIyear, HWlstyear, HWstartmembyear, HWendmembyear, ndayexedthresholdyear, DDthresholdyear, subHWarrayyear, fitsubHWarrayyear, sstMeanarrayyear, impossible_fit = calc_HWMIyear(data[ilon,ilat,:,:,:], cross_valid = cv, percent_thresh = percent_thresh, duration_min = duration_min)\n",
    "        #print HWMIyear\n",
    "        HWMI[:,:,ilat,ilon]=np.array(HWMIyear)\n",
    "        ndayexedthreshold[:,:,ilat,ilon]=np.array(ndayexedthresholdyear)\n",
    "        subHWarray[:,:, :, ilat,ilon]=np.array(subHWarrayyear)\n",
    "        #DDthreshold[:,:,ilat,ilon]=np.array(DDthresholdyear)\n",
    "        #fitsubHWarray[:,:, :, ilat,ilon]=np.array(fitsubHWarrayyear)\n",
    "        #sstMeanarray[:,:,ilat,ilon]=np.array(sstMeanarrayyear)\n",
    "        #impossible_fit_list.append((impossible_fit,ilat,ilon))\n",
    "        return()\n",
    "\n",
    "Parallel(n_jobs=-1, timeout = 5*3600, verbose = 20, require='sharedmem', mmap_mode='w+')(delayed(parallelized_HWMIs_computation)(ilat, ilon) for ilat in range(nlat) for ilon in range(nlon))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "built-small",
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.isdir(dirout)==False:\n",
    "    os.makedirs(dirout)\n",
    "years=range(ref_year1,ref_year2)\n",
    "j=0\n",
    "for i,iyear in enumerate(range(ref_year1, ref_year2+1)):\n",
    "    #print iyear\n",
    "    parameters_str = reg_name+\"_\"+season+\"_\"+cv_str+'_percent%i'%(percent_thresh)+'_daymin%i'%(duration_min)+\"_ref_\"+str(ref_year1)+\"_\"+str(ref_year2)+\"_year_\"+str(iyear)\n",
    "    varout1 = \"HWMI\"+\"_\"+var+\"_\"+parameters_str\n",
    "    vout1=\"HWMI\"+\"_\"+var\n",
    "    fileout=dirout+varout1+\".nc\"#0%i01, monstart)\n",
    "    print(fileout)\n",
    "    if len(glob(fileout))==1:\n",
    "        os.remove(fileout)\n",
    "    fout=netCDF4.Dataset(fileout, \"w\")\n",
    "    #fin=netCDF4.Dataset(targetflst[iyear])\n",
    "    lat = fout.createDimension('lat', nlat)\n",
    "    lon = fout.createDimension('lon', nlon) \n",
    "    rea = fout.createDimension('realisation', nrealisation) \n",
    "    latitudes = fout.createVariable('lat', np.float32, ('lat',)) \n",
    "    longitudes = fout.createVariable('lon', np.float32,  ('lon',)) \n",
    "\n",
    "    # Time variable\n",
    "    timedim = fout.createDimension('time', None) \n",
    "    times = fout.createVariable('time', np.float64, ('time',)) \n",
    "    times.units = \"hours since 1970-01-01 00:00:00\"\n",
    "    times.calendar = 'gregorian' \n",
    "    times[:]=date2num(datetime(iyear,season_start_day[0],season_start_day[1]), units=times.units,calendar = times.calendar) ####______HERE_____\n",
    "\n",
    "\n",
    "    latitudes[:] = lats_reg\n",
    "    lonaux = lons_reg\n",
    "    longitudes[:] = lonaux\n",
    "    latitudes.units = 'degree_north'  \n",
    "    longitudes.units = 'degree_east' \n",
    "\n",
    "    # Create the HWMI 4-d variable\n",
    "    HWMIfile = fout.createVariable(vout1, np.float32, ('time','realisation','lat','lon'))\n",
    "    fout.description = 'HWMI index (Russo et al. 2014) for ' + season + ' computed in cross validation'\n",
    "    fout.history = 'computed from python script by C.Prodhomme & S.Lecestre' + time.ctime(time.time())\n",
    "    fout.source = 'HWMI for ' + expname\n",
    "    latitudes.units = 'degree_north'  \n",
    "    longitudes.units = 'degree_east' \n",
    "    HWMIfile.units = 'Probability'\n",
    "\n",
    "    # Create the nb of days 4-d variable\n",
    "    expercentfile = fout.createVariable(\"nbdaygtpercentpct_\"+var, np.float32, ('time','realisation','lat','lon')) \n",
    "    expercentfile.units = 'Number of days'  \n",
    "\n",
    "    # Write the HWMI variable\n",
    "    HWMIaux=HWMI[i:i+1,0:1,:,:]\n",
    "    HWMIfile[0:1,0:1,:,:]=HWMIaux\n",
    "\n",
    "    # Write the number of days\n",
    "    exedaux=ndayexedthreshold[i:i+1,j:j+1,:,:]\n",
    "    expercentfile[0:1,0:1,:,:]=exedaux \n",
    "\n",
    "    # Write the DDthreshold\n",
    "    #DDthresholdfile=DDthreshold\n",
    "    fout.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e02c2cf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}