From 8dbdb8dac484cbc300f98a20646b89a25f06820b Mon Sep 17 00:00:00 2001 From: simLC <84915517+simLC@users.noreply.github.com> Date: Fri, 4 Jun 2021 15:14:16 +0200 Subject: [PATCH] Upload of HWMI and HWs_detection packages HWMI.py was already needed for Calc-HWMI.ipynb script. HWs_detection.py and Detect_HWs_script.ipynb working versions. --- Detect_HWs_script.ipynb | 312 ++++++++++++++++++++++++++++++++++ HWMI.py | 250 +++++++++++++++++++++++++++ HWs_detection.py | 365 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 927 insertions(+) create mode 100644 Detect_HWs_script.ipynb create mode 100644 HWMI.py create mode 100644 HWs_detection.py diff --git a/Detect_HWs_script.ipynb b/Detect_HWs_script.ipynb new file mode 100644 index 0000000..038ecd4 --- /dev/null +++ b/Detect_HWs_script.ipynb @@ -0,0 +1,312 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib\n", + "#matplotlib.use('Agg')" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import datetime\n", + "import time\n", + "import copy\n", + "import shutil\n", + "import sys\n", + "sys.path.append('/cnrm/pastel/USERS/lecestres/analyse/')\n", + "from function_read import *\n", + "from HWs_detection import *\n", + "import numpy as np\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", + "import matplotlib.pyplot as plt\n", + "from mpl_toolkits.basemap import Basemap, shiftgrid\n", + "\n", + "import calendar\n", + "import locale\n", + "locale.setlocale( locale.LC_ALL , 'en_US' )\n", + "from netCDF4 import num2date, date2num \n", + "import netCDF4\n", + "from math import sin, cos, sqrt, atan2, radians\n", + "from sklearn.neighbors import DistanceMetric\n", + "from math import radians\n", + "\n", + "from joblib import Parallel, delayed\n", + "import joblib" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "### EXPERIENCE NAME\n", + "expname = \"ocean_reanalysis_GREP\"\n", + "#expname = \"sst_retroprevision_sys7\"\n", + "\n", + "### PERCENTILE THRESHOLD\n", + "percent_thresh = 95\n", + "#percent_thresh = 90\n", + "\n", + "### MINIMAL DURATION OF A HW\n", + "duration_min = 5\n", + "#duration_min = 3\n", + "\n", + "### REGION OF EXPERIENCE\n", + "reg_name = \"global\"\n", + "#reg_name = 'lat30_70-lon-40_0'\n", + "#reg_name = 'lat-20_30-lon-30_10'\n", + "if reg_name == 'global':\n", + " lats_bnds = [-90,90]\n", + " lons_bnds = [-180,180] \n", + "elif reg_name == 'lat30_70-lon-40_0':\n", + " lats_bnds = [30,70]\n", + " lons_bnds = [-40, 0]\n", + "elif reg_name == 'lat-20_30-lon-30_10':\n", + " lats_bnds = [-20, 20]\n", + " lons_bnds = [-30, 10]\n", + "nlat = lats_bnds[1]-lats_bnds[0]\n", + "nlon = lons_bnds[1]-lons_bnds[0]\n", + "\n", + "### YEARS\n", + "if expname == 'ocean_reanalysis_GREP':\n", + " end_year=2018\n", + " start_year=1993\n", + "elif expname == 'sst_retroprevision_sys7':\n", + " end_year=2016\n", + " start_year=1993\n", + "nyear=end_year-start_year+1\n", + "\n", + "### SEASON\n", + "season = \"NDJFMAM\"\n", + "if season == 'NDJFMAM':\n", + " nday = 211\n", + " season_start_day = [11,1] #1stNov\n", + " season_end_day = [5,31] #31stMay\n", + "elif season == 'DJF':\n", + " nday = 90\n", + "ndayseas = nday//duration_min +1\n", + "\n", + "if expname == \"ocean_reanalysis_GREP\": \n", + " ### NUMBER OF MEMBS\n", + " nmemb = 1\n", + " chosen_memb = [0] #should never change\n", + " memb_str = 'memb' + str(chosen_memb[0])\n", + " \n", + " ### CROSS VALIDATION\n", + " cv = False\n", + " if cv:\n", + " cv_str = \"CV\"\n", + " else:\n", + " cv_str = 'notCV'\n", + " \n", + "elif expname == \"sst_retroprevision_sys7\": \n", + " ### NUMBER OF MEMBS\n", + " nmemb = 1 #For now keep just one memb by one\n", + " chosen_memb = [21] #to be filled in numerical order\n", + " if len(chosen_memb) != nmemb:\n", + " raise(ValueError)\n", + " if nmemb>1:\n", + " memb_str = 'memb' + str(chosen_memb[0]) + '-' + str(chosen_memb[-1])\n", + " else:\n", + " memb_str = 'memb' + str(chosen_memb[0])\n", + " \n", + " ### CROSS VALIDATION\n", + " cv = True\n", + " if cv:\n", + " cv_str = \"CV\"\n", + " else:\n", + " cv_str = 'notCV'\n", + "\n", + "### PARAMETERS\n", + "parameters_str = reg_name + \"_\" + season + \"_\" + cv_str + '_percent%i'%(percent_thresh) + '_daymin%i'%(duration_min) + \"ref%i-%i\"%(start_year, end_year) + '_' + memb_str\n", + "\n", + "### NAME OF THE VARIABLES IN THE NC FILES\n", + "#varname = 'HWMI_' + parameters_str\n", + "varname = 'subHW_' + parameters_str" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "pathHWMI = '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/' + expname + '/' + memb_str + '/' + varname + '/'\n", + "files = glob(pathHWMI + '*.nc')\n", + "files.sort()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1993_1994.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1994_1995.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1995_1996.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1996_1997.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1997_1998.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1998_1999.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1999_2000.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2000_2001.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2001_2002.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2002_2003.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2003_2004.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2004_2005.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2005_2006.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2006_2007.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2007_2008.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2008_2009.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2009_2010.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2010_2011.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2011_2012.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2012_2013.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2013_2014.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2014_2015.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2015_2016.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2016_2017.nc',\n", + " '/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/HWMI_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_2017_2018.nc']" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "files" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n", + "varf.variables[varname].shape : (1, 1, 180, 360)\n" + ] + } + ], + "source": [ + "obsregyearslst = []\n", + "for file in files:\n", + " varf=netCDF4.Dataset(file)\n", + " varf.variables[varname]\n", + " vararray, lats_reg, lons_reg = extract_array(varf, varname, ndayseas, np.array(lons_bnds), np.array(lats_bnds), start_time = 0, level=0)\n", + " obsreg = vararray[:, np.newaxis, :, :]\n", + " maskobs = obsreg.mask\n", + " obsregyearslst.append(obsreg)\n", + "\n", + "obsregyears = np.ma.array(obsregyearslst)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n", + "(array([0, 0, 0, ..., 0, 0, 0]), array([ 13, 13, 13, ..., 175, 175, 175]), array([ 0, 1, 2, ..., 218, 219, 220]))\n", + "/cnrm/pastel/USERS/lecestres/NO_SAVE/data/ocean_reanalysis_GREP/memb0/Ampli_Fields_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0/Ampli_Fields_global_NDJFMAM_notCV_percent95_daymin5ref1993-2018_memb0_1993_1994.nc\n", + "ndayseas : 1\n", + "(434, 1)\n", + "(434, 1)\n", + "time for year 1993 : 3.8596889972686768\n", + "(array([0, 0, 0, ..., 0, 0, 0]), array([ 13, 13, 13, ..., 174, 174, 174]), array([ 0, 1, 2, ..., 230, 231, 232]))\n" + ] + }, + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mstart_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mexpname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreg_name\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmemb_str\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparameters_str\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstart_year\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlats_reg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlons_reg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mHWampliobsmembyear\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mHWstartobsmembyear\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mHWendobsmembyear\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfieldobslstmembyear\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalc_HW_MY\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobsregyears\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaskobs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlats_reg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlons_reg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mallowdist\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Total time for detection : '\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mstart_time\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/mnt/nfs/d50/pastel/USERS/lecestres/analyse/HWs_detection.py\u001b[0m in \u001b[0;36mcalc_HW_MY\u001b[0;34m(mod, mask, lat, lon, args, allowdist)\u001b[0m\n\u001b[1;32m 271\u001b[0m \u001b[0;31m#print(lat)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 272\u001b[0m \u001b[0;31m#print(lon)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 273\u001b[0;31m \u001b[0mHW\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdetectHW1year\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfield\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlon\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mallowdist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 274\u001b[0m \u001b[0;31m#print('HW list :', HW)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 275\u001b[0m \u001b[0mstart_time_i\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/mnt/nfs/d50/pastel/USERS/lecestres/analyse/HWs_detection.py\u001b[0m in \u001b[0;36mdetectHW1year\u001b[0;34m(field, lat, lon, args, allowdist)\u001b[0m\n\u001b[1;32m 84\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mnei\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlistnei\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 85\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnei\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mptHWlst\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m#Not interested if neighbour has already been looked for\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 86\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mnei\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mHWpointaux\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m#if neighbour point is indeed part of the HW\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 87\u001b[0m \u001b[0;31m#add the neighbourg to the seedlist\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 88\u001b[0m \u001b[0mseedlist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnei\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], + "source": [ + "start_time = time.time()\n", + "args = (expname, reg_name, memb_str, parameters_str, start_year, lats_reg, lons_reg)\n", + "HWampliobsmembyear, HWstartobsmembyear, HWendobsmembyear, fieldobslstmembyear = calc_HW_MY(obsregyears, maskobs, lats_reg, lons_reg, args, allowdist=1)\n", + "print('Total time for detection : ', time.time() - start_time)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/HWMI.py b/HWMI.py new file mode 100644 index 0000000..e571b58 --- /dev/null +++ b/HWMI.py @@ -0,0 +1,250 @@ +import numpy as np +from function_read import * +from scipy.stats import gaussian_kde + +def runningmean(field, ndrun): + nday = field.shape[0] + run=np.zeros((nday)) + run.shape + run = np.convolve(field, np.ones((ndrun,))/ndrun)[(ndrun-1-ndrun/2):(-ndrun+ndrun/2)] + return(run) + + +def find_subHW(target, percent_thresh, duration_min, cross_valid, opt="polyfit"): + """function that find all the sub heatwave from an array + target: numpy array of 3dimension: nyear, nday, members + percent_thresh: Chosen percentile (usually 90, 95 or 99) + duration_min: Minimum length of a HW or a sub HW + return: all sub heatwaves in a list of list of list and ecfmax to construct the fit""" + + #Do we have interest in keeping cross validation? Yes we have + + targetaux=np.array(target) + #print(target.shape) + (nyear, nday, nmemb) = target.shape + + if cross_valid: + #if cross_validation, only use of the first memb + perc_thresh_climlst=[] + for iyear in range(nyear): + #print(targetaux.shape) + #targetaux=np.delete(targetaux[:,:,0], iyear, axis=0) + perc_thresh_clim_not_smoothed=np.percentile(np.delete(targetaux[:,:,0], iyear, axis=0), percent_thresh, axis=0) + if opt=="runningmean": + ndrun=30 + perc_thresh_clim=runningmean(perc_thresh_clim_not_smoothed,ndrun) + #compute the smoothed climatology using a polyfit of order 6 + if opt=="polyfit": + perc_thresh_clim = np.polyval(np.polyfit(range(nday), perc_thresh_clim_not_smoothed, deg=6), range(nday)) + perc_thresh_climlst.append(perc_thresh_clim) + perc_thresh_clim=np.array(perc_thresh_climlst) + perc_thresh_clim=extend_table(perc_thresh_clim, [nmemb]) + perc_thresh_clim=np.swapaxes(np.swapaxes(perc_thresh_clim,1,0), 2,1) + + else: + targetaux2 = np.swapaxes(targetaux, 1,2) + targetaux2 = np.reshape(targetaux2, (nyear * nmemb, nday)) + perc_thresh_clim_not_smoothed = np.percentile(targetaux2[:,:], percent_thresh, axis=0) + if opt=="runningmean": + ndrun=30 + perc_thresh_clim = runningmean(perc_thresh_clim_not_smoothed,ndrun) + #compute the smoothed climatology using a polyfit of order 6 + if opt=="polyfit": + perc_thresh_clim = np.polyval(np.polyfit(range(nday), perc_thresh_clim_not_smoothed, deg=6), range(nday)) + perc_thresh_clim = extend_table(perc_thresh_clim, [nyear, nmemb]) + perc_thresh_clim = np.swapaxes(perc_thresh_clim, 1,2) + + #array to keep subHW values + ndayseas = nday//duration_min + 1 + subHWarray = np.zeros((nyear, ndayseas, nmemb)) + + #check days exceeding the percentile + #exed_percentile = (target[:,limseas1:limseas2, :] - perc_thresh_clim[:, limseas1:limseas2, :]) + exed_percentile = (target - perc_thresh_clim) + exed_percentile = np.ma.array(exed_percentile, mask=exed_percentile<0) + + ndaysexed_percentile = np.sum(exed_percentile>0, axis = 1) + DD_threshold = np.ma.sum(exed_percentile, axis = 1) #np.ma.sum?, np.ma? + + #find heatwave (>x days exceeding threshold percentile) + mask_duration_min = np.array(1-exed_percentile.mask, copy=True) + #print('mask_duration_min.shape : ', mask_duration_min.shape) + for ilen in range(1,duration_min): + mask_duration_min[:,0:-(ilen),:] = mask_duration_min[:,0:-(ilen),:] + (1-exed_percentile.mask[:,ilen:,:]) + #remove heat wave shorter than duration_min + mask_duration_min[mask_duration_min 0: + subHWlst.append(subHW) + subHWarray[iyear, itime//duration_min+(ilen-1)//duration_min+1, imemb] = subHW + #print("itime",itime) + + #keep the subHWlst list in the HW list + HWlst.append(subHWlst) + # keep the beginning and end of the HW + HWstart.append(itime) + HWend.append(itime + ilen + duration_min-1) + #keep the maximum to fit the ecdf + maxtmp.append(max(subHWlst)) + #keep the max for the considered season + if len(maxtmp) != 0: + if imemb == 0: + ecdfmax.append(max(maxtmp)) + HWlstmemb.append(HWlst) + HWstartmemb.append(HWstart) + HWendmemb.append(HWend) + else: + #print imemb, ilat, ilon + #print "no heatwave" + if imemb==0: + ecdfmax.append(-1) + HWlstmemb.append([[0]]) + HWstartmemb.append([0]) + HWendmemb.append([0]) + #keep HW for years + HWlstmembyear.append(HWlstmemb) + HWstartmembyear.append(HWstartmemb) + HWendmembyear.append(HWendmemb) + + return(HWlstmembyear, HWstartmembyear, HWendmembyear, ecdfmax, ndaysexed_percentile, DD_threshold, subHWarray) + + +def fitforHWMI(ecdfmax, cross_valid, yeartoexclude = None): + """compute the gaussian fit + ecdfmax: list of yearly maximum of heatwave for function find_subHW + cross_valid: True/False + year_to_exclude: int (only if using cross validation + """ + ecdfmax=np.array(ecdfmax) + if cross_valid: + ecdfmax = np.delete(ecdfmax, yeartoexclude) + + #remove -1 from ecdf max (corresponding to years with no heatwave) + ecdfmaxfit = ecdfmax[ecdfmax>-1] + try: + kde = gaussian_kde(ecdfmaxfit) + except: + #print("impossible to fit")# : lon %i, lat%i"%(ilon,ilat)) + #print ecdfmax + kde = False + return(kde) + +def calcHWMImemb(HWlstmemb, subHWarrayyear, nday, duration_min, nmemb, kde): + """ + calc the HWMI for a given year + HWlstmemb: list of list containing the list of subheatwave for each members + kde: gaussian fit from the function fitforHWMI + return: list """ + #print HWlstmemb + HWMI = [] + HWmembfit = [] + fitsubHWarray = np.zeros((nday//duration_min+1, nmemb)) + #print nmemb + for imemb in range(nmemb): + HWlstfit = [] + maxHWnorm=0 + HWlst=HWlstmemb[imemb] + #print HWlst + #compute HW magnitude based on normalized sub heatwave magnitude + HWlstfit = [] + #print mask_duration_min[iyear,:,imemb,ilat,ilon] + for iHW in range(len(HWlst)): + sublst=HWlst[iHW] + #find the problability corresponding to the subHW on the fitted cdf + #print sublst + if kde!=False: + fitsum=sum([kde.integrate_box_1d(0,i) for i in sublst]) + #print fit + #keep the maximum + maxHWnorm=max(fitsum,maxHWnorm) + #keep all the heatwave + HWlstfit.append(fitsum) + #feet the subHW from the array to keep all subHW values (not optimal but easier) + for imemb in range(nmemb): + for iday in range(nday//duration_min+1): + subaux = subHWarrayyear[iday, imemb] + if subaux>0: + fitsubHWarray[iday, imemb] = kde.integrate_box_1d(0,subaux) + #print imemb, iHW, fitsum + #store HW magnitude + HWMI.append(maxHWnorm) + HWmembfit.append(HWlstfit) + #print HWMI + return(HWMI, HWmembfit, fitsubHWarray) + +def calc_HWMIyear(target, cross_valid, percent_thresh, duration_min): + """ + function computing the HWMI index for a 3D array at a certain_loc ilat, ilon + target: numpy array of 3dimension: nyear, nday, nmembs + cross_valid: True/False, in case of True, exclude the year + return the HWMI for each year + """ + + (nyear,nday,nmemb) = target.shape + HWlstmembyear, HWstartmembyear, HWendmembyear, ecdfmax, ndaysexed_percentile, DD_threshold, subHWarrayyear = find_subHW(target, percent_thresh=percent_thresh, duration_min=duration_min, cross_valid=cross_valid) + HWMIyear = [] + HWlstyear = [] + fitsubHWarrayyear = np.zeros((nyear, nday//duration_min+1, nmemb)) + count_impossible_fit = 0 + if not(cross_valid): + kde = fitforHWMI(ecdfmax, cross_valid = False) + if kde == False: + count_impossible_fit = 1 + for iyear in range(nyear): + if cross_valid: + kde = fitforHWMI(ecdfmax, cross_valid = True, yeartoexclude = iyear) + #print HWlstmembyear[iyear] + #print len(HWlstmembyear[iyear]) + #exit(1) + #print(subHWarrayyear.shape) + if kde == False: + count_impossible_fit += 1/nyear + HWMI, HWlstfit, fitsubHWarray = calcHWMImemb(HWlstmembyear[iyear], subHWarrayyear[iyear, :, :], nday, duration_min, nmemb, kde) + HWMIyear.append(HWMI) + HWlstyear.append(HWlstfit) + fitsubHWarrayyear[iyear, :, :] = fitsubHWarray + + return(HWMIyear, HWlstyear, HWstartmembyear, HWendmembyear, ndaysexed_percentile, DD_threshold, subHWarrayyear, fitsubHWarrayyear, count_impossible_fit) \ No newline at end of file diff --git a/HWs_detection.py b/HWs_detection.py new file mode 100644 index 0000000..72cae52 --- /dev/null +++ b/HWs_detection.py @@ -0,0 +1,365 @@ +import numpy as np +import os +import datetime +import time +import copy +import shutil +import sys +from function_read import * +import joblib +from joblib import Parallel, delayed + +def detectHW1year(field, lat, lon, args, allowdist=1): + """ + field : np.array fitsubHW + lat, lon : np.arrays corresponding to lat, lon range. + allowdist : neighbourhood geometrical radius. The temporal radius is fixed to one. + """ + expname, reg_name, memb_str, parameters_str, start_year, lats_reg, lons_reg = args + #print("lon",lon) + #print(lat) + nlat= len(lat) + nlon=len(lon) + #print(field.shape) + #print(field) + HWwhere = np.ma.where(field>0) #select indices ix of field where field[ix]>0 --> There is a HW at ix + #Maybe field>0.05 would be better? + + print(HWwhere) + nHWpoints = HWwhere[0].shape[0] #number of points with a HW + + #transform HWwhere in a list of points + HWpoint = [] + for iHW in range(nHWpoints): + HWpoint.append((HWwhere[0][iHW],HWwhere[1][iHW], HWwhere[2][iHW])) + # 0 --> time variable : day + # 1,2 --> space variable : lat, lon + + # + #_______sort heatwave points by neigbours________ + # + + HWpointaux=list(HWpoint) #make a copy + HW = [] + iHW = 0 + iyear=0 + #initialize the list of seeds with the first point + seedlist = [HWpointaux[0]] + #remove seed from the list + HWpointaux=HWpointaux[1:] + #print seedlist + + #run over all the points + while len(HWpointaux)>0: #still some points we did not reach + + #create a list to store the points of one HW + ptHWlst = [] + while len(seedlist)>0: + #print(seedlist) + #remove the seed from the list of seeds and keep the current seed point + seedpoint=seedlist[0] + #print seedlist + #add the seed to the heatwave + #print ptHWlst + #check neighbours for spatial and temporal dimensions + listnei = spacelist_neighbors_highdist2(seedpoint, allowdist) + + # adding temporal neighbours + neibef = (seedpoint[0]-1, seedpoint[1], seedpoint[2]) + + #if not(neibef in listnei): #&(dist>0): + # listnei.append(neibef) + + neiaft = (seedpoint[0]+1, seedpoint[1], seedpoint[2]) + #if not(neiaft in listnei): #&(dist>0): + # listnei.append(neiaft) + + listnei = listnei+[neibef, neiaft] + + if reg_name != "global": + #remove element outside the limits (avoid useless parcours of HWpointaux) + listnei = [nei for nei in listnei if all(0 <= x < y for x, y in zip(nei, field.shape))] + #need to have lats_range between 0 and 179 et lons_range between 0 and 359 + + for nei in listnei: + if not(nei in ptHWlst): #Not interested if neighbour has already been looked for + if nei in HWpointaux: #if neighbour point is indeed part of the HW + #add the neighbourg to the seedlist + seedlist.append(nei) + #remove the neigbourg from the heatwave list + HWpointaux.remove(nei) + #print(seedlist) + #add point to HW list + ptHWlst.append(seedpoint) + # + seedlist=seedlist[1:] + + #once the seed list is empty take a new seed it remains points + #print HWpointaux + if len(HWpointaux)>0: + seedlist = [HWpointaux[0]] + HWpointaux=HWpointaux[1:] + #keep the list of point for each HW + HW.append(ptHWlst) + + # + #_______END of sort heatwave points by neigbours________ + # + + return HW + +def computeHWcara(field, mask, lat, lon, HW): + """ + field : np.array of shape (ndayseas, nlat, nlon) + """ + #print mask.shape + # + #create containers to stored Heatwave informations + HWampli = [] + HWstart = [] + HWend = [] + fieldHWlst = [] + #mask3d = np.repeat(mask[np.newaxis,:,:], field.shape[0], axis=0) + #print(mask3d.shape) + #mask3d.shape = field.shape + #loop over all the HW + #for pointlst in HW: + for iHW in range(len(HW)): + #print("Dealing with new heatwave") + pointlst = HW[iHW] + #loop over all the point part of the HW + # unmask HW point only over land + start = 10000 + end = -1 + #create a fully masked field + fieldHW = np.zeros(field.shape) + for point in pointlst: + #exclude oceanic points for HW calculations + #print(iHW) + #print(point) + #print(mask.shape) + #if not(mask[point[1], point[2]]): + fieldHW[point[0], point[1], point[2]] = field[point[0], point[1], point[2]] + + # keep heatwave values only if some land points have been found + #print(fieldHW.shape, mask.shape) + fieldHW = np.ma.array(fieldHW, mask=mask[:,0, :,:]) #) + #fieldHW = np.ma.array(fieldHW) #) + ampli = area_av(fieldHW, 1, 2, lat, lon, opt="mean") + + #print(np.ma.where(fieldHW>0)) + if np.ma.where(fieldHW>0)[0].size != 0: + start = np.min(np.ma.where(fieldHW>0)[0]) + end = np.max(np.ma.where(fieldHW>0)[0]) + else: + start = None + end = None + #print('ampli : ', ampli) + #check if HW is only oceanic + + if np.sum(ampli)>0: + HWampli.append(ampli) + HWstart.append(start) #start date + HWend.append(end) #end date + fieldHWlst.append(np.ma.sum(fieldHW, axis=0)) + #else: + #print("Ce test n est pas verifie") + return(HWampli, HWstart, HWend, fieldHWlst) + +def calc_HW_MY(mod, mask, lat, lon, args, allowdist=1): + """ + mod : data np.array of shape (nyear, ndayseas, nmemb, nlon, nlat) + """ + ###_______DEF INSIDE FUNCTION + + def fonction_export_allHWs(iyear, nHW_i, dirout, HWampliobs_yeari, fieldobslst_yeari, args): + nlat, nlon = fieldobslst_yeari[0].shape + expname, reg_name, memb_str, parameters_str, start_year, lats_reg, lons_reg = args + fileout= dirout+'Ampli_Fields_'+ parameters_str +"_%i_%i.nc"%(iyear+start_year,iyear+start_year+1) + if len(glob(fileout))==1: + os.remove(fileout) + print(fileout) + fout=netCDF4.Dataset(fileout, "w") + + # Create Dimensions + lat = fout.createDimension('lat', nlat) + lon = fout.createDimension('lon', nlon) + #rea = fout.createDimension('realization', nmemb) + timedim = fout.createDimension('time', None) + HWdim = fout.createDimension('nHW', nHW_i) #To adapt for each year + days = fout.createDimension('days', ndayseas) # Check with chloe, maybe use time instead + # Create Variables + times = fout.createVariable('time', np.float64, ('time',)) + latitudes = fout.createVariable('lat', np.float32, ('lat',)) + longitudes = fout.createVariable('lon', np.float32, ('lon',)) + HWvar = fout.createVariable('nHW', np.int_, ('nHW',)) + daysvar = fout.createVariable('days', np.float32, ('days')) + + # Time variable + # Useless, keep just in case + times.units = 'hours since 0001-01-01 00:00:00' + times.calendar = 'gregorian' + times[:]=date2num(datetime((start_year+iyear),12,1), units = times.units, calendar = times.calendar) + + # Fill Variables + latitudes[:] = lats_reg + lonaux = lons_reg + lonaux[lonaux<0]=lonaux[lonaux<0]+360 + longitudes[:] = lonaux + latitudes.units = 'degree_north' + longitudes.units = 'degree_east' + + # WHAT VARIABLE DO WE WANT TO EXPORT? + # Create Export Variables + + + Amplifile = fout.createVariable('Ampli', np.float32, ('nHW', 'days')) + Fieldfile = fout.createVariable('Field', np.float32, ('nHW', 'lat', 'lon')) + + fout.description = 'HW Amplitude and Fields computed from HWMI indexes files' #find something better with Chloe + fout.history = 'computed from python script by C.Prodhomme and S.Lecestre' + time.ctime(time.time()) + fout.source = 'HW Amplitude and Fields for '+ expname + latitudes.units = 'degree_north' + longitudes.units = 'degree_east' + + Amplifile.units = 'Amplitude per days' + Fieldfile.units = 'Amplitude per surface' + + + + # Write Ampli and Fields + # To do so, transform list of d-arrays into (d+1)-arrays + #maskFalse = np.zeros((nHW_i, ndayseas), dtype = bool) + #ampliaux = np.ma.array(np.zeros((nHW_i, ndayseas)), mask = maskFalse) + print('ndayseas : ', ndayseas) + ampliaux = np.zeros((nHW_i, ndayseas)) + maskFalse2 = np.zeros((nHW_i, nlat, nlon), dtype = bool) + fieldaux = np.ma.array(np.zeros((nHW_i, nlat, nlon)), mask = maskFalse2) + for iHW in range(nHW_i): + #print(ampliaux[iHW, :].shape, np.array(HWampliobs_yeari[iHW][:]).shape) + #print('HWampliobs_yeari[0].shape : ', HWampliobs_yeari[0].shape) + #print(HWampliobs_yeari[0]) + ampliaux[iHW, :] = np.array(HWampliobs_yeari[iHW][:]) + #ampliaux[iHW, :].mask = (ampliaux[iHW, :]==0) #keep only days with ampli non equal to 0 + fieldaux[iHW, :, :] = np.array(fieldobslst_yeari[iHW][:,:]) + fieldaux[iHW, :, :].mask = np.array(fieldobslst_yeari[iHW][:,:].mask) #np.array to make an indpt copy + print(Amplifile[:,:].shape) + print(ampliaux[:,:].shape) + Amplifile[:,:] = ampliaux[:,:] + Fieldfile[:,:,:] = fieldaux[:,:,:] + + fout.close() + return + + ### END OF INSIDE FUNCTION + + nyear, ndayseas, nmemb, nlon, nlat = mod.shape + expname, reg_name, memb_str, parameters_str, start_year, lats_reg, lons_reg = args + HWamplimembyear = [] + HWstartmembyear = [] + HWendmembyear = [] + fieldlstmembyear = [] + for imemb in range(nmemb): + print(imemb) + HWampliyear = [] + HWstartyear = [] + HWendyear = [] + fieldlstyear = [] + for iyear in range(nyear): + #print(iyear) + field=mod[iyear,:, imemb,:,:] + #print(lat) + #print(lon) + HW = detectHW1year(field, lat, lon, args, allowdist) + #print('HW list :', HW) + start_time_i = time.time() + HWampli, HWstart, HWend, fieldHWlst = computeHWcara(field, mask, lat, lon, HW) + HWampliyear.append(HWampli) + HWstartyear.append(HWstart) + HWendyear.append(HWend) + fieldlstyear.append(fieldHWlst) + #c'est ici qu'on vient rajouter des bails pour un bail à un seul membre + #fontion_plot(iyear, dirout) + nHW_i = len(HW) + dirout = "/cnrm/pastel/USERS/lecestres/NO_SAVE/data/"+ expname +'/' + memb_str + '/Ampli_Fields_'+ parameters_str +'/' + if not(os.path.isdir(dirout)): + os.makedirs(dirout) + fonction_export_allHWs(iyear, nHW_i, dirout, HWampli, fieldHWlst, args) + print('time for year ', start_year + iyear,' : ', time.time()-start_time_i) + HWamplimembyear.append(HWampliyear) + HWstartmembyear.append(HWstartyear) + HWendmembyear.append(HWendyear) + fieldlstmembyear.append(fieldlstyear) + + + return(HWamplimembyear, HWstartmembyear, HWendmembyear, fieldlstmembyear) + +def spacelist_neighbors_highdist(point, allowdist=3): + """ + point: list of length 3 : (time,lat,lon) + """ + if allowdist == 1: + neisouth = (point[0], (point[1]-1), (point[2] +180)%360 - 180) + neinorth = (point[0], (point[1]+1), (point[2] +180)%360 - 180) + neiwest = (point[0], (point[1]), (point[2]-1 +180)%360 - 180) + neieast = (point[0], (point[1]), (point[2]+1 +180)%360 - 180) + neiNW = (point[0], (point[1]+1), (point[2]-1 +180)%360 - 180) + neiNE = (point[0], (point[1]+1), (point[2]+1 +180)%360 - 180) + neiSW = (point[0], (point[1]-1), (point[2]-1 +180)%360 - 180) + neiSE = (point[0], (point[1]-1), (point[2]+1 +180)%360 - 180) + listnei = [neisouth, neinorth, neiwest, neieast, neiNW, neiNE, neiSW, neiSE] + + else: + spaceneighbours = [] + allowdistx=allowdist+1 + allowdisty=allowdist+1 + for idistx in range(allowdistx): + for idisty in range(allowdisty): + # adding space neighbours + nei1 = (point[0], (point[1]-idisty), (point[2]-idistx + 180)%360 - 180) + nei2 = (point[0], (point[1]-idisty), (point[2]+idistx + 180)%360 - 180) + nei3 = (point[0], (point[1]+idisty), (point[2]+idistx + 180)%360 - 180) + nei4 = (point[0], (point[1]+idisty), (point[2]-idistx + 180)%360 - 180) + spaceneighbours += [nei1, nei2, nei3, nei4] + listnei = list(set(listnei)) #sort and leave duplicates terms + if point in listnei: + #should be always true + listnei.remove(point) + + return(listnei) + + + +def spacelist_neighbors_highdist2(point, allowdist=3): + """ + point: list of length 3 : (time,lat,lon) + """ + if allowdist == 1: + neisouth = (point[0], (point[1]-1), point[2]%360) + neinorth = (point[0], (point[1]+1), point[2]%360) + neiwest = (point[0], (point[1]), (point[2]-1)%360) + neieast = (point[0], (point[1]), (point[2]+1)%360) + neiNW = (point[0], (point[1]+1), (point[2]-1)%360) + neiNE = (point[0], (point[1]+1), (point[2]+1)%360) + neiSW = (point[0], (point[1]-1), (point[2]-1)%360) + neiSE = (point[0], (point[1]-1), (point[2]+1)%360) + listnei = [neisouth, neinorth, neiwest, neieast, neiNW, neiNE, neiSW, neiSE] + + else: + spaceneighbours = [] + allowdistx=allowdist+1 + allowdisty=allowdist+1 + for idistx in range(allowdistx): + for idisty in range(allowdisty): + # adding space neighbours + nei1 = (point[0], (point[1]-idisty), (point[2]-idistx)%360) + nei2 = (point[0], (point[1]-idisty), (point[2]+idistx)%360) + nei3 = (point[0], (point[1]+idisty), (point[2]+idistx)%360) + nei4 = (point[0], (point[1]+idisty), (point[2]-idistx)%360) + spaceneighbours += [nei1, nei2, nei3, nei4] + listnei = list(set(listnei)) #sort and leave duplicates terms + if point in listnei: + #should be always true + listnei.remove(point) + + return(listnei) \ No newline at end of file