Experiments/GHCN/sphere_GHCN_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DAILY GLOBAL HISTORICAL CLIMATOLOGY NETWORK (GHCN-DAILY) \n",
    "\n",
    "Version 3.24\n",
    "\n",
    "[Andreas Loukas](https://andreasloukas.wordpress.com/), [EPFL LTS2](https://lts2.epfl.ch/), [Michaël Defferrard](http://deff.ch)\n",
    "\n",
    "Menne, M.J., I. Durre, B. Korzeniewski, S. McNeal, K. Thomas, X. Yin, S. Anthony, R. Ray, \n",
    "R.S. Vose, B.E.Gleason, and T.G. Houston, 2012: Global Historical Climatology Network - \n",
    "Daily (GHCN-Daily), Version 3. \n",
    "\n",
    "The five core features are:\n",
    "\n",
    "* PRCP = Precipitation (tenths of mm)\n",
    "* SNOW = Snowfall (mm)\n",
    "* SNWD = Snow depth (mm)\n",
    "* TMAX = Maximum temperature (tenths of degrees C)\n",
    "* TMIN = Minimum temperature (tenths of degrees C)\n",
    "\n",
    "The other features are:\n",
    "\n",
    "* ACMC = Average cloudiness midnight to midnight from 30-second ceilometer data (percent)\n",
    "* ACMH = Average cloudiness midnight to midnight from manual observations (percent)\n",
    "* ACSC = Average cloudiness sunrise to sunset from 30-second ceilometer data (percent)\n",
    "* ACSH = Average cloudiness sunrise to sunset from manual observations (percent)\n",
    "* AWDR = Average daily wind direction (degrees)\n",
    "* AWND = Average daily wind speed (tenths of meters per second)\n",
    "* DAEV = Number of days included in the multiday evaporation total (MDEV)\n",
    "* DAPR = Number of days included in the multiday precipiation total (MDPR)\n",
    "* DASF = Number of days included in the multiday snowfall total (MDSF)\t\t  \n",
    "* DATN = Number of days included in the multiday minimum temperature (MDTN)\n",
    "* DATX = Number of days included in the multiday maximum temperature (MDTX)\n",
    "* DAWM = Number of days included in the multiday wind movement (MDWM)\n",
    "* DWPR = Number of days with non-zero precipitation included in multiday precipitation total (MDPR)\n",
    "* EVAP = Evaporation of water from evaporation pan (tenths of mm)\n",
    "* FMTM = Time of fastest mile or fastest 1-minute wind (hours and minutes, i.e., HHMM)\n",
    "* FRGB = Base of frozen ground layer (cm)\n",
    "* FRGT = Top of frozen ground layer (cm)\n",
    "* FRTH = Thickness of frozen ground layer (cm)\n",
    "* GAHT = Difference between river and gauge height (cm)\n",
    "* MDEV = Multiday evaporation total (tenths of mm; use with DAEV)\n",
    "* MDPR = Multiday precipitation total (tenths of mm; use with DAPR and DWPR, if available)\n",
    "* MDSF = Multiday snowfall total \n",
    "* MDTN = Multiday minimum temperature (tenths of degrees C; use with DATN)\n",
    "* MDTX = Multiday maximum temperature (tenths of degress C; use with DATX)\n",
    "* MDWM = Multiday wind movement (km)\n",
    "* MNPN = Daily minimum temperature of water in an evaporation pan (tenths of degrees C)\n",
    "* MXPN = Daily maximum temperature of water in an evaporation pan (tenths of degrees C)\n",
    "* PGTM = Peak gust time (hours and minutes, i.e., HHMM)\n",
    "* PSUN = Daily percent of possible sunshine (percent)\n",
    "* SN*# = Minimum soil temperature (tenths of degrees C) where * corresponds to a code for ground cover and # corresponds to a code for soil depth. For ground cover and depth codes see [readme](https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt)\n",
    "* SX*# = Maximum soil temperature (tenths of degrees C)  where * corresponds to a code for ground cover  and # corresponds to a code for soil depth.  See SN*# for ground cover and depth codes. \n",
    "* TAVG = Average temperature (tenths of degrees C) [Note that TAVG from source 'S' corresponds to an average for the period ending at 2400 UTC rather than local midnight]\n",
    "* THIC = Thickness of ice on water (tenths of mm)\t\n",
    "* TOBS = Temperature at the time of observation (tenths of degrees C)\n",
    "* TSUN = Daily total sunshine (minutes)\n",
    "* WDF1 = Direction of fastest 1-minute wind (degrees)\n",
    "* WDF2 = Direction of fastest 2-minute wind (degrees)\n",
    "* WDF5 = Direction of fastest 5-second wind (degrees)\n",
    "* WDFG = Direction of peak wind gust (degrees)\n",
    "* WDFI = Direction of highest instantaneous wind (degrees)\n",
    "* WDFM = Fastest mile wind direction (degrees)\n",
    "* WDMV = 24-hour wind movement (km)\t   \n",
    "* WESD = Water equivalent of snow on the ground (tenths of mm)\n",
    "* WESF = Water equivalent of snowfall (tenths of mm)\n",
    "* WSF1 = Fastest 1-minute wind speed (tenths of meters per second)\n",
    "* WSF2 = Fastest 2-minute wind speed (tenths of meters per second)\n",
    "* WSF5 = Fastest 5-second wind speed (tenths of meters per second)\n",
    "* WSFG = Peak gust wind speed (tenths of meters per second)\n",
    "* WSFI = Highest instantaneous wind speed (tenths of meters per second)\n",
    "* WSFM = Fastest mile wind speed (tenths of meters per second)\n",
    "* WT** = Weather Type. For numeric codes ** see [readme](https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt)\n",
    "* WV** = Weather in the Vicinity. For numeric codes ** see [readme](https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt)\n",
    "\n",
    "[FTP link](ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import sys\n",
    "sys.path.append('../..')\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"1\"  # change to chosen GPU to use, nothing if work on CPU\n",
    "\n",
    "import numpy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import healpy as hp\n",
    "import pandas as pd\n",
    "\n",
    "from tqdm import tqdm\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deepsphere import models, experiment_helper, plot, utils\n",
    "from deepsphere.data import LabeledDatasetWithNoise, LabeledDataset\n",
    "import hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = \"/mnt/nas/LTS2/datasets/ghcn-daily/processed/\"\n",
    "rawpath = \"/mnt/nas/LTS2/datasets/ghcn-daily/raw/\"\n",
    "newdatapath = \"./data/ghcn-daily/processed/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "years = np.arange(2010,2015)\n",
    "years2 = np.arange(2010,2018)\n",
    "\n",
    "feature_names = ['PRCP', 'TMIN', 'TMAX', 'SNOW']\n",
    "n_features = len(feature_names)\n",
    "n_years  = len(years)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Station"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'stations_{:4d}-{:4d}.npz'.format(years[0], years[-1])\n",
    "\n",
    "# only recompute if necessary\n",
    "if not os.path.isfile(datapath+filename):\n",
    "    \n",
    "    # Variable   Columns   Type\n",
    "    # ID            1-11   Character\n",
    "    # LATITUDE     13-20   Real\n",
    "    # LONGITUDE    22-30   Real\n",
    "    # ELEVATION    32-37   Real\n",
    "    # STATE        39-40   Character\n",
    "    # NAME         42-71   Character\n",
    "    # GSN FLAG     73-75   Character\n",
    "    # HCN/CRN FLAG 77-79   Character\n",
    "    # WMO ID       81-85   Character\n",
    "\n",
    "    id_ghcn, lat, lon, elev, name = [], [], [], [], []\n",
    "    with open(rawpath+'ghcnd-stations.txt', 'r') as f:\n",
    "        for line in f:\n",
    "\n",
    "            iid, ilat, ilon, ielev, iname = line[0:11], line[12:20], line[21:30], line[31:37], line[41:71]         \n",
    "            # state, GSN_flag, HCNCRN_flag, WMO_id = line[38:40], line[72:75], line[76:79], line[80:85] \n",
    "\n",
    "            assert (not iid.isspace()) and (not ilat.isspace()) and (not ilon.isspace()) \\\n",
    "                and (not ielev.isspace()) and (not iname.isspace())\n",
    "\n",
    "            id_ghcn.append(iid.strip())\n",
    "            lat.append(float(ilat.strip()))\n",
    "            lon.append(float(ilon.strip()))\n",
    "            elev.append(float(ielev.strip()))\n",
    "            name.append(iname.strip())\n",
    "\n",
    "    id_ghcn, lat, lon, elev, name = np.array(id_ghcn), np.array(lat), np.array(lon), np.array(elev), np.array(name)\n",
    "\n",
    "    # Identify *relevant* stations: These are stations giving measurements in the years of interest\n",
    "    # that also have known coordinates.\n",
    "\n",
    "    # first, construct the set of all ghcn identifiers encountered in daily\n",
    "    id_ghcn_relevant = set([])\n",
    "\n",
    "    for yearIdx,year in enumerate(years):\n",
    "\n",
    "        filename2 = rawpath+'{:4}.csv.gz'.format(year)\n",
    "        print('- pre-parsing : {}'.format(filename2))\n",
    "\n",
    "        df = pd.read_csv(filename2, names=['id_ghcn', 'date', 'type', 'value', '?0', '?1', '?2', '?3'], \\\n",
    "                         nrows=None, usecols=[0,1,2,3])\n",
    "\n",
    "\n",
    "        id_ghcn_relevant |= set(df[\"id_ghcn\"].values)\n",
    "\n",
    "    # second, find identifiers both in id_ghcn and id_ghcn_relevant\n",
    "    id_ghcn_relevant = set(id_ghcn) & id_ghcn_relevant\n",
    "\n",
    "    # third, keep only relevant station data \n",
    "    keep = [id in id_ghcn_relevant for id in id_ghcn] \n",
    "    id_ghcn, lat, lon, elev, name = id_ghcn[keep], lat[keep], lon[keep], elev[keep], name[keep] \n",
    "\n",
    "    # free up some memory\n",
    "    del id_ghcn_relevant, keep\n",
    "\n",
    "    np.savez_compressed(newdatapath+filename, id_ghcn=id_ghcn, lat=lat, lon=lon, elev=elev, name=name, years=years)\n",
    "    \n",
    "else:\n",
    "    station_file = np.load(datapath+filename)\n",
    "    id_ghcn, lat, lon, elev, name = station_file['id_ghcn'], station_file['lat'], station_file['lon'], station_file['elev'], station_file['name']\n",
    "    del station_file\n",
    "    \n",
    "n_stations = id_ghcn.shape[0]\n",
    "print('{} weather stations identified.'.format(n_stations))\n",
    "\n",
    "#  a dictionary mapping GHCN ids to local ids (rows in id array) \n",
    "ghcn_to_local = dict(zip(id_ghcn, np.arange(n_stations)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = ['PRCP', 'TMIN', 'TMAX', 'SNOW', 'SNWD', 'AWND', 'WT']\n",
    "# max stations, 50'469\n",
    "# in 2010, in nbr stations [36k, 16k, 16k, 6k, 22k, 20k, 0, 0, 28, 1k, 6, 5k, 0]\n",
    "# in 2011, in nbr stations [36k, 16k, 16k, 6k, 22k, 18k, 0, 0, 38, 1k, 22, 5k, 0]\n",
    "# in 2014, in nbr stations [39k, 15k, 15k, 6k, 24k, 19k, 0, 0, 40, 1k, 49, 4k, 0]\n",
    "\n",
    "filenames = []\n",
    "datas = []\n",
    "valid_days_list = []\n",
    "n_stations_list = []\n",
    "for feature_name in feature_names:\n",
    "    filenames.append('data_{:4d}-{:4d}_{}.npz'.format(years[0], years[-1], feature_name))\n",
    "    print(f'- Checking if file {filenames[-1]} exists..')\n",
    "\n",
    "    # only recompute if necessary\n",
    "    if not os.path.isfile(newdatapath+filenames[-1]):\n",
    "\n",
    "        print('- The file is not there. Parsing everything from raw. This will take a while.')\n",
    "        os.makedirs(newdatapath, exist_ok=True)\n",
    "        # Load the station measurements into a year-list of dataframes\n",
    "        df_years = []\n",
    "\n",
    "        for yearIdx,year in enumerate(years):\n",
    "\n",
    "            filename_year = rawpath+'{:4}.csv.gz'.format(year)\n",
    "            print(' - parsing *{}*'.format(filename_year))\n",
    "\n",
    "            df = pd.read_csv(filename_year, names=['id_ghcn', 'date', 'type', 'value', 'MF', 'qualityF', 'source', '?0'], \\\n",
    "                             nrows=None, usecols=[0,1,2,3,5])\n",
    "\n",
    "            # create a new column with the id_local\n",
    "            id_local = [ghcn_to_local.get(id_g) for id_g in df[\"id_ghcn\"].values]\n",
    "            id_local = [-1 if v is None else v for v in id_local]\n",
    "            id_local = np.array(id_local).astype(np.int)\n",
    "\n",
    "            df = df.assign(id_local=pd.Series(id_local, index=df.index).values)\n",
    "\n",
    "            # remove measurement of stations with unknown id_local\n",
    "            df = df[df.id_local != -1] \n",
    "            \n",
    "            # replace measurements with bad quality flag\n",
    "            #df.value[~df.qualityF.isna()] = np.nan\n",
    "            df = df[df.qualityF.isna()]\n",
    "            df = df.drop('qualityF', axis=1)\n",
    "            \n",
    "            df_years.append(df)\n",
    "\n",
    "        del df, id_local\n",
    "        print('done!')\n",
    "\n",
    "        # Construct one array per feature and save it to disk\n",
    "\n",
    "        # indicate for which days we have measurements (this also deals with months of different lengths)\n",
    "        valid_days = np.zeros((n_years, 12, 31), dtype=np.bool)\n",
    "\n",
    "        for _, name in enumerate(feature_names):\n",
    "\n",
    "            print(f' - Looking at {name}')\n",
    "\n",
    "            data = np.zeros((n_stations, n_years, 12, 31), dtype=np.float) * np.nan\n",
    "\n",
    "            for yearIdx,year in enumerate(years):\n",
    "\n",
    "                df = df_years[yearIdx]\n",
    "                idf = df.loc[df.type.str.contains(name)]\n",
    "\n",
    "                print(f'  - year {year}')\n",
    "\n",
    "                # remove measurement of stations with unknown id_local\n",
    "                idf = idf[idf.id_local != -1] \n",
    "\n",
    "                for monthIdx,month in enumerate(range(1,12+1)): \n",
    "                    for dayIdx,day in enumerate(range(1,31+1)):        \n",
    "\n",
    "                        date = int('{:4d}{:02d}{:02d}'.format(year,month,day))\n",
    "                        jdf = idf.loc[idf['date'] == date]\n",
    "\n",
    "                        # sort data according to the id_local \n",
    "                        jdf.set_index('id_local', inplace=True)\n",
    "                        jdf = jdf.sort_index()\n",
    "\n",
    "                        index = jdf.index.values\n",
    "                        if name is 'WT' or name is 'WV':\n",
    "                            values = jdf.type.str.extract(r'(\\d+)').values.astype(int)\n",
    "                            values = values[:,0]\n",
    "                        else:\n",
    "                            values = jdf['value'].values.astype(np.float)\n",
    "\n",
    "                        if len(index) != 0: \n",
    "                            data[index,yearIdx,monthIdx,dayIdx] = values\n",
    "                            valid_days[yearIdx,monthIdx,dayIdx] = True\n",
    "\n",
    "            print('  - saving to disk')\n",
    "            np.savez_compressed(newdatapath+'data_{:4d}-{:4d}_{}.npz'.format(years[0], years[-1], name), data=data, valid_days=valid_days)\n",
    "\n",
    "            del index, values, df, idf, jdf    \n",
    "\n",
    "    else:    \n",
    "        print('- Loading data from disk..')\n",
    "\n",
    "        data_file = np.load(newdatapath+filenames[-1])\n",
    "        data, valid_days = data_file['data'], data_file['valid_days']        \n",
    "        n_stations = data.shape[0]\n",
    "        print(f'- {n_stations} stations loaded.')\n",
    "        data = data.reshape((n_stations, n_years*12*31))\n",
    "        if feature_name == 'TMIN' or feature_name == 'TMAX' or feature_name == 'PRCP':\n",
    "            data = data.astype(np.float)\n",
    "            data /= 10\n",
    "        datas.append(data)\n",
    "        valid_days = np.squeeze(valid_days.reshape(n_years*12*31)).astype(np.bool)\n",
    "        valid_days_list.append(valid_days)\n",
    "        n_stations_list.append(n_stations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#feature_name = 'TMAX'\n",
    "# feature_names\n",
    "filenames = []\n",
    "datas_old = []\n",
    "valid_days_list = []\n",
    "n_stations_list = []\n",
    "for feature_name in feature_names[:4]:\n",
    "    filenames.append('data_{:4d}-{:4d}_{}.npz'.format(years[0], years[-1], feature_name))\n",
    "    print(f'- Checking if file {filenames[-1]} exists..')\n",
    "\n",
    "    # only recompute if necessary\n",
    "    if not os.path.isfile(datapath+filenames[-1]):\n",
    "\n",
    "        print('- The file is not there. Parsing everything from raw. This will take a while.')\n",
    "        # Load the station measurements into a year-list of dataframes\n",
    "        df_years = []\n",
    "\n",
    "        for yearIdx,year in enumerate(years):\n",
    "\n",
    "            filename_year = rawpath+'{:4}.csv.gz'.format(year)\n",
    "            print(' - parsing *{}*'.format(filename_year))\n",
    "\n",
    "            df = pd.read_csv(filename_year, names=['id_ghcn', 'date', 'type', 'value', '?0', '?1', '?2', '?3'], \\\n",
    "                             nrows=None, usecols=[0,1,2,3])\n",
    "\n",
    "            # create a new column with the id_local\n",
    "            id_local = [ghcn_to_local.get(id_g) for id_g in df[\"id_ghcn\"].values]\n",
    "            id_local = [-1 if v is None else v for v in id_local]\n",
    "            id_local = np.array(id_local).astype(np.int)\n",
    "\n",
    "            df = df.assign(id_local=pd.Series(id_local, index=df.index).values)\n",
    "\n",
    "            # remove measurement of stations with unknown id_local\n",
    "            df = df[df.id_local != -1] \n",
    "\n",
    "            df_years.append(df)\n",
    "\n",
    "        del df, id_local\n",
    "        print('done!')\n",
    "\n",
    "        # Construct one array per feature and save it to disk\n",
    "\n",
    "        # indicate for which days we have measurements (this also deals with months of different lengths)\n",
    "        valid_days = np.zeros((n_years, 12, 31), dtype=np.bool)\n",
    "\n",
    "        for _, name in enumerate(feature_names):\n",
    "\n",
    "            print(f' - Looking at {name}')\n",
    "\n",
    "            data = np.zeros((n_stations, n_years, 12, 31), dtype=np.float) * np.nan\n",
    "\n",
    "            for yearIdx,year in enumerate(years):\n",
    "\n",
    "                df = df_years[yearIdx]\n",
    "                idf = df.loc[df['type'] == name]\n",
    "\n",
    "                print(f'  - year {year}')\n",
    "\n",
    "                # remove measurement of stations with unknown id_local\n",
    "                idf = idf[idf.id_local != -1] \n",
    "\n",
    "                for monthIdx,month in enumerate(range(1,12+1)): \n",
    "                    for dayIdx,day in enumerate(range(1,31+1)):        \n",
    "\n",
    "                        date = int('{:4d}{:02d}{:02d}'.format(year,month,day))\n",
    "                        jdf = idf.loc[idf['date'] == date]\n",
    "\n",
    "                        # sort data according to the id_local \n",
    "                        jdf.set_index('id_local', inplace=True)\n",
    "                        jdf = jdf.sort_index()\n",
    "\n",
    "                        index = jdf.index.values\n",
    "                        values = jdf['value'].values.astype(np.float)\n",
    "\n",
    "                        if len(index) != 0: \n",
    "                            data[index,yearIdx,monthIdx,dayIdx] = values\n",
    "                            valid_days[yearIdx,monthIdx,dayIdx] = True\n",
    "\n",
    "            print('  - saving to disk')\n",
    "            np.savez_compressed('processed/data_{:4d}-{:4d}_{}.npz'.format(years[0], years[-1], name), data=data, valid_days=valid_days)\n",
    "\n",
    "            del index, values, df, idf, jdf    \n",
    "\n",
    "    else:    \n",
    "        print('- Loading data from disk..')\n",
    "\n",
    "        data_file = np.load(datapath+filenames[-1])\n",
    "        data, valid_days = data_file['data'], data_file['valid_days']        \n",
    "        n_stations = data.shape[0]\n",
    "        print(f'- {n_stations} stations loaded.')\n",
    "        data = data.reshape((n_stations, n_years*12*31))\n",
    "        if feature_name == 'TMIN' or feature_name == 'TMAX' or feature_name == 'PRCP':\n",
    "            data = data.astype(np.float)\n",
    "            data /= 10\n",
    "        datas_old.append(data)\n",
    "        valid_days = np.squeeze(valid_days.reshape(n_years*12*31)).astype(np.bool)\n",
    "        valid_days_list.append(valid_days)\n",
    "        n_stations_list.append(n_stations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(full_data[:,:,1], full_data[:,:,2], s=2)\n",
    "plt.plot(range(-40, 41), range(-40, 41), 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(datas[1][:,valid_days], datas[2][:,valid_days])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert n_stations_list[0] == n_stations_list[1] == n_stations_list[2] == n_stations_list[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.all(valid_days_list[0] == valid_days_list[0])\n",
    "assert np.all(valid_days_list[0] == valid_days_list[1])\n",
    "assert np.all(valid_days_list[0] == valid_days_list[2])\n",
    "assert np.all(valid_days_list[0] == valid_days_list[3])\n",
    "# assert np.all(valid_days_list[0] == valid_days_list[4])\n",
    "# assert np.all(valid_days_list[0] == valid_days_list[5])\n",
    "# assert np.all(valid_days_list[0] == valid_days_list[6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_data = np.stack(datas_old, axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_data = full_data[:, valid_days_list[0], :]\n",
    "\n",
    "n_days = full_data.shape[1]\n",
    "\n",
    "assert n_stations == full_data.shape[0]\n",
    "\n",
    "print(f'n_stations: {n_stations}, n_days: {n_days}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.all(np.isnan(datas[0]).all(axis=(1)) == np.isnan(datas[1]).all(axis=(1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keep = ~np.isnan(full_data).all(axis=1) # keep a station if we have at least one measurement in the period ot interest \n",
    "keep = np.all(keep, axis=1)\n",
    "\n",
    "data = full_data[keep,:, :]\n",
    "n_stations, n_days, _ = data.shape\n",
    "print(f'n_stations: {n_stations}, n_days: {n_days}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keep_days = ~np.isnan(data).all(axis=0)\n",
    "keep_days = np.all(keep_days, axis=1)\n",
    "print(keep_days.shape)\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_data = full_data[:,:,1:3]\n",
    "keep_temp = ~np.isnan(temp_data).all(axis=1)\n",
    "keep_temp = np.all(keep_temp, axis=1)\n",
    "\n",
    "temp_data = temp_data[keep_temp,:, :]\n",
    "n_stations, n_days, _ = temp_data.shape\n",
    "print(f'n_stations: {n_stations}, n_days: {n_days}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n",
    "\n",
    "ax.set_global()\n",
    "# ax.stock_img()\n",
    "ax.coastlines()\n",
    "\n",
    "plt.plot(lon[keep], lat[keep], 'or', marker='o', markerfacecolor='r', markersize=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "\n",
    "size = 40\n",
    "fig = plt.figure(figsize=(2*size, 1*size))\n",
    "for i in range(len(feature_names)):\n",
    "    ax = fig.add_subplot(4, 1, i+1, projection=ccrs.PlateCarree())\n",
    "    ax.set_global()\n",
    "    ax.coastlines()\n",
    "    ax.scatter(lon[keep], lat[keep], s=None, c=data[:,0,i], cmap=plt.get_cmap('RdYlGn'))\n",
    "    ax.set_title(feature_names[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seems to have outliers data. Can be sure to be correct (Tmax = -99 when Tmin = -8) or (Tmin = -72.8 for multiple Tmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = np.where(data<-70)\n",
    "#print(data[out[:2]])\n",
    "data[out] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[:,:,2], data[:,:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[:,:,2], data[:,:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[9,:,2], data[9,:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[:,:,1], data[:,:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[:,:,2], data[:,:,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[:,:,0], data[:,:,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(data[3,:,0], data[3,:,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 10))\n",
    "plt.plot(data[:6,:,2].T, 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 10))\n",
    "plt.plot(data[0,:,0], 'o')\n",
    "plt.figure(figsize=(20, 10))\n",
    "plt.plot(data[0,:,2], 'o')\n",
    "plt.figure(figsize=(20, 10))\n",
    "plt.plot(data[1,:,2], 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 10))\n",
    "plt.plot(data[3,:,0], 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pygsp import utils\n",
    "from pygsp.graphs import Graph\n",
    "from pygsp.graphs import NNGraph\n",
    "class sphereGraph(NNGraph):\n",
    "    def __init__(self, phi, theta, neighbors, rad=True, epsilon=True, **kwargs):\n",
    "        if not rad:\n",
    "            theta, phi = np.deg2rad(theta), np.deg2rad(phi)\n",
    "        theta -= np.pi/2\n",
    "        ct = np.cos(theta).flatten()\n",
    "        st = np.sin(theta).flatten()\n",
    "        cp = np.cos(phi).flatten()\n",
    "        sp = np.sin(phi).flatten()\n",
    "        x = st * cp\n",
    "        y = st * sp\n",
    "        z = ct\n",
    "        self.coords = np.vstack([x, y, z]).T\n",
    "        NNtype = 'radius' if epsilon else 'knn'\n",
    "        plotting = {\"limits\": np.array([-1, 1, -1, 1, -1, 1])*0.5}\n",
    "        self.n_vertices = len(self.coords)\n",
    "        super(sphereGraph, self).__init__(self.coords, k=neighbors, NNtype=NNtype, center=False, rescale=False,\n",
    "                                     plotting=plotting, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sphereGraph(lon[keep], lat[keep], 100, rad=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(25,25))\n",
    "axes = fig.add_subplot(111, projection='3d')\n",
    "g.plot(vertex_size=10, edges=False, ax=axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.compute_laplacian(\"normalized\")\n",
    "g.compute_fourier_basis(recompute=True, n_eigenvectors=1000)\n",
    "g.set_coordinates(g.U[:,1:4])\n",
    "g.plot(vertex_size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(g.e[:16], 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_full = sphereGraph(lon, lat, 100, rad=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(25,25))\n",
    "axes = fig.add_subplot(111, projection='3d')\n",
    "g_full.plot(vertex_size=10, edges=False, ax=axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_full.compute_laplacian(\"normalized\")\n",
    "g_full.compute_fourier_basis(recompute=True, n_eigenvectors=500)\n",
    "g_full.set_coordinates(g_full.U[:,1:4])\n",
    "g_full.plot(vertex_size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(g_full.e[:16], 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## inference train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "params = {'L': [g.L.astype(np.float32)]*4,\n",
    "          'p': [1,1,1,1],\n",
    "          'F': [10, 20, 50, 1],\n",
    "          'K': [5]*4,\n",
    "          'batch_norm': [True]*4}\n",
    "params['dir_name'] = 'GHCN_essai'\n",
    "params['num_feat_in'] = 1\n",
    "params['conv'] = 'chebyshev5'\n",
    "params['pool'] = 'max'\n",
    "params['activation'] = 'relu'\n",
    "params['statistics'] = None#'mean'\n",
    "params['regularization'] = 0\n",
    "params['dropout'] = 1\n",
    "params['num_epochs'] = 50  # Number of passes through the training data.\n",
    "params['batch_size'] = 32\n",
    "params['scheduler'] = lambda step: tf.train.exponential_decay(5e-1, step, decay_steps=5, decay_rate=1)\n",
    "params['optimizer'] = lambda lr: tf.train.GradientDescentOptimizer(lr)\n",
    "n_evaluations = 200\n",
    "params['eval_frequency'] = int(params['num_epochs'] * n_days / params['batch_size'] / n_evaluations)\n",
    "params['M'] = []\n",
    "params['regression']=True\n",
    "model = models.cgcnn(**params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from deepsphere.data import LabeledDatasetWithNoise, LabeledDataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Temperature MAX from Temperature MIN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = full_data.transpose((1, 0, 2))\n",
    "keepToo = ~np.isnan(dataset[:,:,1:3]).any(axis=0)\n",
    "keepToo = keepToo.all(axis=1)\n",
    "dataset_temp = dataset[:, keepToo, 1:3]\n",
    "\n",
    "# remove outliers (remove outiler stations) (other option is to replace value by mean of knn)\n",
    "out = np.where((dataset_temp<-60) + (dataset_temp>80))\n",
    "dataset_temp = np.delete(dataset_temp, np.unique(out[1]), axis=1)\n",
    "keepToo[np.argwhere(keepToo)[np.unique(out[1])]] = False\n",
    "\n",
    "out = np.where(dataset_temp[:,:,1]==60)\n",
    "dataset_temp = np.delete(dataset_temp, np.unique(out[1]), axis=1)\n",
    "keepToo[np.argwhere(keepToo)[np.unique(out[1])]] = False\n",
    "\n",
    "out = np.where((dataset_temp[:,:,0]==-35.6)*(dataset_temp[:,:,1]>0))\n",
    "dataset_temp = np.delete(dataset_temp, np.unique(out[1]), axis=1)\n",
    "keepToo[np.argwhere(keepToo)[np.unique(out[1])]] = False\n",
    "\n",
    "# replace by means of neighbours (use graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dataset[:,keepToo,1:3].shape)\n",
    "print(np.isnan(dataset[:,keepToo,1:3]).all())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n",
    "\n",
    "ax.set_global()\n",
    "# ax.stock_img()\n",
    "ax.coastlines()\n",
    "\n",
    "plt.plot(lon[keepToo], lat[keepToo], 'or', marker='o', markerfacecolor='r', markersize=2)\n",
    "\n",
    "fig2 = plt.figure(figsize=(25,25))\n",
    "axes = fig2.add_subplot(111, projection='3d')\n",
    "gTemp = sphereGraph(lon[keepToo], lat[keepToo], 10, rad=False)\n",
    "gTemp.plot(vertex_size=10, edges=True, ax=axes)\n",
    "gTemp.compute_laplacian(\"normalized\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## temp max from temp min\n",
    "\n",
    "ratio = 0.7\n",
    "n_days = dataset_temp.shape[0]\n",
    "limit= int(ratio*n_days)\n",
    "\n",
    "mean = dataset_temp.mean(axis=(0,1))[0]\n",
    "std = dataset_temp.std(axis=(0,1))[0]\n",
    "\n",
    "x_train = (dataset_temp[:limit,:,0] - mean) / std\n",
    "labels_train = dataset_temp[:limit,:,1]\n",
    "x_val = (dataset_temp[limit:,:,0] - mean) / std\n",
    "labels_val = dataset_temp[limit:,:,1]\n",
    "\n",
    "\n",
    "training = LabeledDataset(x_train, labels_train)\n",
    "validation = LabeledDataset(x_val, labels_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_val, labels_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train, labels_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Precipirtation from temperatures (MIN, MAX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = full_data.transpose((1, 0, 2)).copy()\n",
    "# remove outliers\n",
    "### min temp is impossible\n",
    "out = np.where(dataset[:,:,1]>100)\n",
    "dataset[:,:,1][out] = np.nan\n",
    "\n",
    "### min temp is always the same and max temp doesn't make sens\n",
    "out = np.where(dataset[:,:,1]<-60)\n",
    "dataset[:,:,1][out] = np.nan\n",
    "\n",
    "### max temp is too big ofr the location and the min temp\n",
    "out = np.where(dataset[:,:,2]>60)\n",
    "dataset[:,:,2][out] = np.nan\n",
    "\n",
    "\n",
    "## min is bigger than max\n",
    "out = np.where(dataset[:,:,1]>dataset[:,:,2])\n",
    "dataset[:,:,1:3][out] = np.nan\n",
    "\n",
    "## dif between two days is bigger than 20° (what max is possible)\n",
    "out = np.where(np.abs(np.diff(dataset[:,:,1:3], axis=0))>20)\n",
    "dataset[:,:,1:3][out] = np.nan\n",
    "dataset[:,:,1:3][out[0]+1,out[1], out[2]] = np.nan\n",
    "\n",
    "\n",
    "keepToo = ~np.isnan(dataset[:,:,:3]).any(axis=0)\n",
    "keepSuper = ((~np.isnan(dataset[:,:,:3])).sum(axis=0)>0.75*dataset.shape[0])\n",
    "keepToo = keepToo.all(axis=1)\n",
    "keepSuper = keepSuper.all(axis=1)\n",
    "dataset_prec = dataset[:, keepToo, :3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plop = np.where(dataset[:,:,0]>3000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset[:,:,0][plop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isnan(dataset_prec).any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"n stations minset = {}\".format(keepToo.sum()))\n",
    "print(\"n stations superset = {}\".format(keepSuper.sum()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n",
    "\n",
    "ax.set_global()\n",
    "# ax.stock_img()\n",
    "ax.coastlines()\n",
    "\n",
    "plt.plot(lon[keepToo], lat[keepToo], 'or', marker='o', markerfacecolor='r', markersize=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbour = 1\n",
    "fig = plt.figure(figsize=(25,25))\n",
    "axes = fig.add_subplot(111, projection='3d')\n",
    "gPrec = sphereGraph(lon[keepToo], lat[keepToo], neighbour, rad=False)\n",
    "gPrec.plot(vertex_size=10, edges=False, ax=axes)\n",
    "gPrec.compute_laplacian(\"normalized\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## precipitation from temp\n",
    "\n",
    "ratio = 0.7\n",
    "n_days = dataset_prec.shape[0]\n",
    "limit= int(ratio*n_days)\n",
    "\n",
    "mean = dataset_prec.mean(axis=(0,1))[1:3]\n",
    "std = dataset_prec.std(axis=(0,1))[1:3]\n",
    "# mean_label = dataset_prec.mean(axis=(0,1))[0]\n",
    "# std_label = dataset_prec.std(axis=(0,1))[0]\n",
    "\n",
    "x_train = (dataset_prec[:limit,:,1:3] - mean) / std\n",
    "labels_train = dataset_prec[:limit,:,0]/10\n",
    "x_val = (dataset_prec[limit:,:,1:3] - mean) / std\n",
    "labels_val = dataset_prec[limit:,:,0]/10\n",
    "\n",
    "# location of stations\n",
    "coords_v = np.stack([lon[keepToo], lat[keepToo]], axis=-1)\n",
    "coords_v = (coords_v-coords_v.mean(axis=0))/coords_v.std(axis=0)\n",
    "# altitude of stations\n",
    "alt_v = elev[keepToo]\n",
    "alt_v = (alt_v-alt_v.mean())/alt_v.std()\n",
    "\n",
    "x_train = np.dstack([x_train, np.repeat(coords_v[np.newaxis,:], x_train.shape[0], axis=0),\n",
    "                     np.repeat(alt_v[np.newaxis,:], x_train.shape[0], axis=0),\n",
    "                     np.repeat(w_days[:limit, np.newaxis], x_val.shape[1], axis=1)])\n",
    "x_val = np.dstack([x_val, np.repeat(coords_v[np.newaxis,:], x_val.shape[0], axis=0),\n",
    "                  np.repeat(alt_v[np.newaxis,:], x_val.shape[0], axis=0),\n",
    "                  np.repeat(w_days[limit:, np.newaxis], x_val.shape[1], axis=1)])\n",
    "\n",
    "training = LabeledDataset(x_train, labels_train)\n",
    "validation = LabeledDataset(x_val, labels_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train[:,:,0], labels_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train[:,:,1], labels_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train[:,:,2], labels_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train[:,:,3], labels_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train[:,:,4], labels_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_train[:,:,0], x_train[:,:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = np.empty_like(x_train[:,:,0])\n",
    "time[:,:] = np.arange(time.shape[0])[:,np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(time, labels_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "future temperature to find, regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_days[31+28::12*31]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_months = np.tile(np.repeat(np.arange(12), 31), years[-1]-years[0]+1)[valid_days]\n",
    "w_days = np.tile(np.arange(365),years[-1]-years[0]+1)\n",
    "w_days = np.insert(w_days, (3*365), 365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = full_data.transpose((1, 0, 2)).copy()\n",
    "# remove outliers\n",
    "### min temp is impossible\n",
    "out = np.where(dataset[:,:,1]>100)\n",
    "dataset[:,:,1][out] = np.nan\n",
    "\n",
    "### min temp is always the same and max temp doesn't make sens\n",
    "out = np.where(dataset[:,:,1]<-60)\n",
    "dataset[:,:,1][out] = np.nan\n",
    "\n",
    "### max temp is too big ofr the location and the min temp\n",
    "out = np.where(dataset[:,:,2]>=60)\n",
    "dataset[:,:,2][out] = np.nan\n",
    "\n",
    "## series of exact same values for the same location. seems strange\n",
    "out = np.where((dataset[:,:,1]==-35.6)*(dataset[:,:,2]>0))\n",
    "dataset[:,:,1][out] = np.nan\n",
    "\n",
    "## min is bigger than max\n",
    "out = np.where(dataset[:,:,1]>dataset[:,:,2])\n",
    "dataset[:,:,1:3][out] = np.nan\n",
    "\n",
    "## dif between two days is bigger than 20° (what max is possible)\n",
    "out = np.where(np.abs(np.diff(dataset[:,:,1:3], axis=0))>20)\n",
    "dataset[:,:,1:3][out] = np.nan\n",
    "dataset[:,:,1:3][out[0]+1,out[1], out[2]] = np.nan\n",
    "\n",
    "\n",
    "keepToo = ~np.isnan(dataset[:,:,1:3]).any(axis=0)\n",
    "keepSuper = ((~np.isnan(dataset[:,:,1:3])).sum(axis=0)>0.75*dataset.shape[0])\n",
    "keepToo = keepToo.all(axis=1)\n",
    "keepSuper = keepSuper.all(axis=1)\n",
    "dataset_temp_reg = dataset[:, keepToo, 1:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"n stations minset = {}\".format(keepToo.sum()))\n",
    "print(\"n stations superset = {}\".format(keepSuper.sum()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n",
    "\n",
    "ax.set_global()\n",
    "# ax.stock_img()\n",
    "ax.coastlines()\n",
    "\n",
    "plt.plot(lon[keepToo], lat[keepToo], 'or', marker='o', markerfacecolor='r', markersize=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighbour = 10\n",
    "fig = plt.figure(figsize=(25,25))\n",
    "axes = fig.add_subplot(111, projection='3d')\n",
    "gReg = sphereGraph(lon[keepToo], lat[keepToo], neighbour, rad=False)\n",
    "gReg.plot(vertex_size=10, edges=False, ax=axes)\n",
    "gReg.compute_laplacian(\"normalized\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ~np.isnan(dataset[:,keepSuper,1:3])\n",
    "dataset_super = dataset[:, keepSuper, 1:3] # n_days graphs of superset\n",
    "gRegSuper = sphereGraph(lon[keepSuper], lat[keepSuper], 100, rad=False)\n",
    "gRegSuper.compute_laplacian(\"normalized\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for day in range(dataset_super.shape[0]):\n",
    "    for feat in range(dataset_super.shape[-1]):\n",
    "        dataset_super[day,:,feat] = regression_tikhonov(gRegSuper, dataset_super[day,:,feat], mask[day,:,feat])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isnan(dataset_super).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import sparse\n",
    "\n",
    "def regression_tikhonov(G, y, M, tau=0):\n",
    "    r\"\"\"Solve a regression problem on graph via Tikhonov minimization.\n",
    "    The function solves\n",
    "    .. math:: \\operatorname*{arg min}_x \\| M x - y \\|_2^2 + \\tau \\ x^T L x\n",
    "    if :math:`\\tau > 0`, and\n",
    "    .. math:: \\operatorname*{arg min}_x x^T L x \\ \\text{ s. t. } \\ y = M x\n",
    "    otherwise.\n",
    "    Parameters\n",
    "    ----------\n",
    "    G : :class:`pygsp.graphs.Graph`\n",
    "    y : array, length G.n_vertices\n",
    "        Measurements.\n",
    "    M : array of boolean, length G.n_vertices\n",
    "        Masking vector.\n",
    "    tau : float\n",
    "        Regularization parameter.\n",
    "    Returns\n",
    "    -------\n",
    "    x : array, length G.n_vertices\n",
    "        Recovered values :math:`x`.\n",
    "    Examples\n",
    "    --------\n",
    "    >>> from pygsp import graphs, filters, learning\n",
    "    >>> import matplotlib.pyplot as plt\n",
    "    >>>\n",
    "    >>> G = graphs.Sensor(N=100, seed=42)\n",
    "    >>> G.estimate_lmax()\n",
    "    Create a smooth ground truth signal:\n",
    "    >>> filt = lambda x: 1 / (1 + 10*x)\n",
    "    >>> filt = filters.Filter(G, filt)\n",
    "    >>> rs = np.random.RandomState(42)\n",
    "    >>> signal = filt.analyze(rs.normal(size=G.n_vertices))\n",
    "    Construct a measurement signal from a binary mask:\n",
    "    >>> mask = rs.uniform(0, 1, G.n_vertices) > 0.5\n",
    "    >>> measures = signal.copy()\n",
    "    >>> measures[~mask] = np.nan\n",
    "    Solve the regression problem by reconstructing the signal:\n",
    "    >>> recovery = learning.regression_tikhonov(G, measures, mask, tau=0)\n",
    "    Plot the results:\n",
    "    >>> fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(10, 3))\n",
    "    >>> limits = [signal.min(), signal.max()]\n",
    "    >>> _ = G.plot_signal(signal, ax=ax1, limits=limits, title='Ground truth')\n",
    "    >>> _ = G.plot_signal(measures, ax=ax2, limits=limits, title='Measures')\n",
    "    >>> _ = G.plot_signal(recovery, ax=ax3, limits=limits, title='Recovery')\n",
    "    >>> _ = fig.tight_layout()\n",
    "    \"\"\"\n",
    "\n",
    "    if tau > 0:\n",
    "\n",
    "        y[M == False] = 0\n",
    "\n",
    "        if sparse.issparse(G.L):\n",
    "\n",
    "            def Op(x):\n",
    "                return (M * x.T).T + tau * (G.L.dot(x))\n",
    "\n",
    "            LinearOp = sparse.linalg.LinearOperator([G.N, G.N], Op)\n",
    "\n",
    "            if y.ndim > 1:\n",
    "                sol = np.empty(shape=y.shape)\n",
    "                res = np.empty(shape=y.shape[1])\n",
    "                for i in range(y.shape[1]):\n",
    "                    sol[:, i], res[i] = sparse.linalg.cg(\n",
    "                        LinearOp, y[:, i])\n",
    "            else:\n",
    "                sol, res = sparse.linalg.cg(LinearOp, y)\n",
    "\n",
    "            # TODO: do something with the residual...\n",
    "            return sol\n",
    "\n",
    "        else:\n",
    "\n",
    "            # Creating this matrix may be problematic in term of memory.\n",
    "            # Consider using an operator instead...\n",
    "            if type(G.L).__module__ == np.__name__:\n",
    "                LinearOp = np.diag(M*1) + tau * G.L\n",
    "            return np.linalg.solve(LinearOp, M * y)\n",
    "\n",
    "    else:\n",
    "\n",
    "        if np.prod(M.shape) != G.n_vertices:\n",
    "            raise ValueError(\"M should be of size [G.n_vertices,]\")\n",
    "\n",
    "        indl = M\n",
    "        indu = (M == False)\n",
    "\n",
    "        Luu = G.L[indu, :][:, indu]\n",
    "        Wul = - G.L[indu, :][:, indl]\n",
    "\n",
    "        if sparse.issparse(G.L):\n",
    "            sol_part = sparse.linalg.spsolve(Luu, Wul.dot(y[indl]))\n",
    "        else:\n",
    "            sol_part = np.linalg.solve(Luu, np.matmul(Wul, y[indl]))\n",
    "\n",
    "        sol = y.copy()\n",
    "        sol[indu] = sol_part\n",
    "\n",
    "        return sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords_v = (coords_v-coords_v.mean(axis=0))/coords_v.std(axis=0)\n",
    "plt.plot(coords_v[:,0], coords_v[:,1], 'o', markersize=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords_v[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## temp future (regression)\n",
    "days_pred = 5\n",
    "\n",
    "n_days, n_stations, n_feature= dataset_temp_reg.shape\n",
    "\n",
    "dataset_x = np.vstack([np.roll(dataset_temp_reg, -i, axis=0) for i in range(days_pred)])\n",
    "dataset_x = dataset_x.reshape(days_pred, n_days, n_stations, n_feature).transpose((1,2,3,0))\n",
    "\n",
    "# mean_x = dataset_x.mean(axis=(0,1))\n",
    "# std_x = dataset_x.std(axis=(0,1))\n",
    "\n",
    "# month_x = np.vstack([np.roll(w_months, -i, axis=0) for i in range(days_pred)])\n",
    "# month_x = month_x.reshape(days_pred, n_days).transpose()\n",
    "\n",
    "days_x = np.vstack([np.roll(w_days, -i, axis=0) for i in range(days_pred)])\n",
    "days_x = days_x.reshape(days_pred, n_days).transpose()\n",
    "\n",
    "x_train = dataset_x[:n_days-days_pred,:,0,:]\n",
    "labels_train = dataset_temp_reg[days_pred:,:,0]\n",
    "x_val = dataset_x[:n_days-days_pred,:,1,:]\n",
    "labels_val = dataset_temp_reg[days_pred:,:,1]\n",
    "\n",
    "# x_train = (x_train-mean_x[0])/std_x[0]\n",
    "# x_val = (x_val-mean_x[1])/std_x[1]\n",
    "\n",
    "# location of stations\n",
    "coords_v = np.stack([lon[keepToo], lat[keepToo]], axis=-1)\n",
    "coords_v = (coords_v-coords_v.mean(axis=0))/coords_v.std(axis=0)\n",
    "# altitude of stations\n",
    "alt_v = elev[keepToo]\n",
    "alt_v = (alt_v-alt_v.mean())/alt_v.std()\n",
    "\n",
    "# can include information on which period we are? month?\n",
    "x_train = np.dstack([x_train, \n",
    "#                     np.broadcast_to(month_x[:n_days-days_pred,np.newaxis, :], x_train.shape),\n",
    "                     np.repeat(coords_v[np.newaxis,:], x_train.shape[0],axis=0),\n",
    "                     np.repeat(alt_v[np.newaxis,:], x_train.shape[0],axis=0),\n",
    "                     np.broadcast_to(days_x[:n_days-days_pred,np.newaxis, :], x_train.shape)])\n",
    "x_val = np.dstack([x_val, \n",
    "#                   np.broadcast_to(month_x[:n_days-days_pred,np.newaxis, :], x_val.shape), \n",
    "                   np.repeat(coords_v[np.newaxis,:], x_train.shape[0],axis=0),\n",
    "                   np.repeat(alt_v[np.newaxis,:], x_train.shape[0],axis=0),\n",
    "                   np.broadcast_to(days_x[:n_days-days_pred,np.newaxis, :], x_val.shape)])\n",
    "\n",
    "training = LabeledDataset(x_train, labels_train)\n",
    "validation = LabeledDataset(x_val, labels_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(dataset_temp_reg[:,:,0], dataset_temp_reg[:,:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timereg = np.empty_like(dataset_temp_reg[:,:,0])\n",
    "timereg[:,:] = np.arange(timereg.shape[0])[:,np.newaxis]\n",
    "plt.scatter(timereg, dataset_temp_reg[:,:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "placereg = np.empty_like(dataset_temp_reg[:,:,0])\n",
    "placereg[:,:] = np.arange(placereg.shape[1])\n",
    "plt.scatter(placereg, dataset_temp_reg[:,:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## is it possible to have a difference of 40 degrees between 1 day?  (must fix to +- 20 degrees)\n",
    "plt.scatter(x_train[:,:,days_pred-1], labels_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Misc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_nan = np.isnan(full_data[:,:,0])\n",
    "fig=plt.figure(figsize=(55,100))\n",
    "ax = fig.add_subplot(421)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data[:,:,0]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data[keep_,:,0])\n",
    "ax = fig.add_subplot(422)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "\n",
    "mat_nan = np.isnan(full_data[:,:,1])\n",
    "ax = fig.add_subplot(423)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data[:,:,1]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data[keep_,:,1])\n",
    "ax = fig.add_subplot(424)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "\n",
    "mat_nan = np.isnan(full_data[:,:,2])\n",
    "ax = fig.add_subplot(425)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data[:,:,2]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data[keep_,:,2])\n",
    "ax = fig.add_subplot(426)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "\n",
    "mat_nan = np.isnan(full_data[:,:,3])\n",
    "ax = fig.add_subplot(427)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data[:,:,3]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data[keep_,:,3])\n",
    "ax = fig.add_subplot(428)\n",
    "ax.spy(mat_nan, aspect='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#feature_name = 'TMAX'\n",
    "# feature_names\n",
    "filenames = []\n",
    "datas = []\n",
    "valid_days_list = []\n",
    "n_stations_list = []\n",
    "n_years = len(years2)\n",
    "for feature_name in feature_names:\n",
    "    filenames.append('data_{:4d}-{:4d}_{}.npz'.format(years[0], years2[-1], feature_name))\n",
    "    print(f'- Checking if file {filenames[-1]} exists..')\n",
    "\n",
    "    # only recompute if necessary\n",
    "    if not os.path.isfile(datapath+filenames[-1]):\n",
    "\n",
    "        print('- The file is not there. Parsing everything from raw. This will take a while.')\n",
    "        # Load the station measurements into a year-list of dataframes\n",
    "        df_years = []\n",
    "\n",
    "        for yearIdx,year in enumerate(years):\n",
    "\n",
    "            filename_year = 'raw/{:4}.csv.gz'.format(year)\n",
    "            print(' - parsing *{}*'.format(filename_year))\n",
    "\n",
    "            df = pd.read_csv(filename_year, names=['id_ghcn', 'date', 'type', 'value', '?0', '?1', '?2', '?3'], \\\n",
    "                             nrows=None, usecols=[0,1,2,3])\n",
    "\n",
    "            # create a new column with the id_local\n",
    "            id_local = [ghcn_to_local.get(id_g) for id_g in df[\"id_ghcn\"].values]\n",
    "            id_local = [-1 if v is None else v for v in id_local]\n",
    "            id_local = np.array(id_local).astype(np.int)\n",
    "\n",
    "            df = df.assign(id_local=pd.Series(id_local, index=df.index).values)\n",
    "\n",
    "            # remove measurement of stations with unknown id_local\n",
    "            df = df[df.id_local != -1] \n",
    "\n",
    "            df_years.append(df)\n",
    "\n",
    "        del df, id_local\n",
    "        print('done!')\n",
    "\n",
    "        # Construct one array per feature and save it to disk\n",
    "\n",
    "        # indicate for which days we have measurements (this also deals with months of different lengths)\n",
    "        valid_days = np.zeros((n_years, 12, 31), dtype=np.bool)\n",
    "\n",
    "        for _, name in enumerate(feature_names):\n",
    "\n",
    "            print(f' - Looking at {name}')\n",
    "\n",
    "            data = np.zeros((n_stations, n_years, 12, 31), dtype=np.float) * np.nan\n",
    "\n",
    "            for yearIdx,year in enumerate(years):\n",
    "\n",
    "                df = df_years[yearIdx]\n",
    "                idf = df.loc[df['type'] == name]\n",
    "\n",
    "                print(f'  - year {year}')\n",
    "\n",
    "                # remove measurement of stations with unknown id_local\n",
    "                idf = idf[idf.id_local != -1] \n",
    "\n",
    "                for monthIdx,month in enumerate(range(1,12+1)): \n",
    "                    for dayIdx,day in enumerate(range(1,31+1)):        \n",
    "\n",
    "                        date = int('{:4d}{:02d}{:02d}'.format(year,month,day))\n",
    "                        jdf = idf.loc[idf['date'] == date]\n",
    "\n",
    "                        # sort data according to the id_local \n",
    "                        jdf.set_index('id_local', inplace=True)\n",
    "                        jdf = jdf.sort_index()\n",
    "\n",
    "                        index = jdf.index.values\n",
    "                        values = jdf['value'].values.astype(np.float)\n",
    "\n",
    "                        if len(index) != 0: \n",
    "                            data[index,yearIdx,monthIdx,dayIdx] = values\n",
    "                            valid_days[yearIdx,monthIdx,dayIdx] = True\n",
    "\n",
    "            print('  - saving to disk')\n",
    "            np.savez_compressed('processed/data_{:4d}-{:4d}_{}.npz'.format(years[0], years[-1], name), data=data, valid_days=valid_days)\n",
    "\n",
    "            del index, values, df, idf, jdf    \n",
    "\n",
    "    else:    \n",
    "        print('- Loading data from disk..')\n",
    "\n",
    "        data_file = np.load(datapath+filenames[-1])\n",
    "        data, valid_days = data_file['data'], data_file['valid_days']        \n",
    "        n_stations = data.shape[0]\n",
    "        print(f'- {n_stations} stations loaded.')\n",
    "        data = data.reshape((n_stations, n_years*12*31))\n",
    "        if feature_name == 'TMIN' or feature_name == 'TMAX':\n",
    "            data = data.astype(np.float)\n",
    "            data /= 10\n",
    "        datas.append(data)\n",
    "        valid_days = np.squeeze(valid_days.reshape(n_years*12*31)).astype(np.bool)\n",
    "        valid_days_list.append(valid_days)\n",
    "        n_stations_list.append(n_stations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_data2 = np.stack(datas, axis=2)\n",
    "full_data2 = full_data2[:, valid_days_list[0], :]\n",
    "\n",
    "n_days = full_data2.shape[1]\n",
    "\n",
    "assert n_stations == full_data2.shape[0]\n",
    "\n",
    "print(f'n_stations: {n_stations}, n_days: {n_days}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_nan = np.isnan(full_data2[:,:,0])\n",
    "fig=plt.figure(figsize=(55,100))\n",
    "ax = fig.add_subplot(421)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data2[:,:,0]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data2[keep_,:,0])\n",
    "ax = fig.add_subplot(422)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "\n",
    "mat_nan = np.isnan(full_data2[:,:,1])\n",
    "ax = fig.add_subplot(423)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data2[:,:,1]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data2[keep_,:,1])\n",
    "ax = fig.add_subplot(424)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "\n",
    "mat_nan = np.isnan(full_data2[:,:,2])\n",
    "ax = fig.add_subplot(425)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data2[:,:,2]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data2[keep_,:,2])\n",
    "ax = fig.add_subplot(426)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "\n",
    "mat_nan = np.isnan(full_data2[:,:,3])\n",
    "ax = fig.add_subplot(427)\n",
    "ax.spy(mat_nan, aspect='auto')\n",
    "keep_ = ~np.isnan(full_data2[:,:,3]).all(axis=1)\n",
    "mat_nan = np.isnan(full_data2[keep_,:,3])\n",
    "ax = fig.add_subplot(428)\n",
    "ax.spy(mat_nan, aspect='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sum(~np.isnan(full_data2[:,:,0]).any(axis=1)))\n",
    "print(sum(~np.isnan(full_data2[:,:,1]).any(axis=1)))\n",
    "print(sum(~np.isnan(full_data2[:,:,2]).any(axis=1)))\n",
    "print(sum(~np.isnan(full_data2[:,:,3]).any(axis=1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set1=set(np.where(~np.isnan(full_data2[:,:,0]).any(axis=1))[0])\n",
    "set2=set(np.where(~np.isnan(full_data2[:,:,1]).any(axis=1))[0])\n",
    "set3=set(np.where(~np.isnan(full_data2[:,:,2]).any(axis=1))[0])\n",
    "set4=set(np.where(~np.isnan(full_data2[:,:,3]).any(axis=1))[0])\n",
    "print(len(set1.intersection(set2, set3, set4)))\n",
    "print(len(set2.intersection(set3)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(~np.isnan(full_data[keep_,:,3][:,ploup[0]])).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plop = np.where((~np.isnan(full_data[keep_,:,3][:,ploup[0]])).sum(axis=1)>69)  # stations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ploup = np.where((~np.isnan(full_data[:,:,3])).sum(axis=0)>15000)  # days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_nan = np.isnan(full_data[keep_,:,3][plop[0],:])\n",
    "fig=plt.figure(figsize=(25,25))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.spy(mat_nan, aspect='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(20, 10))\n",
    "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n",
    "\n",
    "ax.set_global()\n",
    "# ax.stock_img()\n",
    "ax.coastlines()\n",
    "\n",
    "plt.plot(lon[plop], lat[plop], 'or', marker='o', markerfacecolor='r', markersize=2)"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}