<jupyter_start><jupyter_text>Climate event detection taskSandbox for preprocessing and first learning test<jupyter_code>%load_ext autoreload
%autoreload 2
%matplotlib inline
import os
import shutil
import sys
sys.path.append('../..')

os.environ["CUDA_VISIBLE_DEVICES"] = "0"  # change to chosen GPU to use, nothing if work on CPU

import numpy as np
import time
import matplotlib.pyplot as plt
import healpy as hp
import pandas as pd

from tqdm import tqdm
from mpl_toolkits.mplot3d import Axes3D
import cartopy.crs as ccrs
import h5py
import matplotlib.patches as patches
from deepsphere import models, experiment_helper, plot, utils
from deepsphere.data import LabeledDatasetWithNoise, LabeledDataset
import hyperparameters<jupyter_output><empty_output><jupyter_text>List of parameters:* PRECT: Total (convective and large-scale) precipitation rate (liq + ice)* PS: Surface pressure* PSL: sea level pressure* QREFHT: Reference height humidity* T200: temp at 200 mbar pressure surface* T500: temp at 500 mbar pressure surface* TMQ: Total (vertically integrated) precipitatable water* TS: Surface temperature (radiative)* U850: Zonal wind at 850 mbar pressure surface* UBOT: Lowest model level zonal wind* V850: Meridional wind at 850 mbar pressure surface* VBOT: Lowest model level meridional wind* Z100: Geopotential Z at 100 mbar pressure surface* Z200: Geopotential Z at 200 mbar pressure surface* ZBOT: Lowest model level height4 measures per day, 365 days a yearresolution of 768 x 1152 equirectangular grid (25-km at equator)boxes:* ymin* xmin* ymax* xmax* class:    * 0: Tropical Depression    * 1: Tropical Cyclone    * 2: Extratropical Cyclone    * 3: Atmospheric River<jupyter_code>year = 2005
datapath = '../../../data/ExtremeWeather/climo_{}.h5'.format(year)
h5f = h5py.File(datapath)
images = h5f["images"] # (1460,16,768,1152) numpy array
boxes = h5f["boxes"] # (1460,15,5) numpy array
lon_ = np.arange(1152)/1152*360
lat_ = np.arange(768)/768*180-90
lon, lat = np.meshgrid(lon_, lat_)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.set_global()
# ax.stock_img()
ax.coastlines()

plt.scatter(lon, lat, s=1,
            c=data[0,:,:], cmap=plt.get_cmap('RdYlBu_r'), alpha=1)

color = ['r', 'y', 'c', 'g']

for i in range(15):
    ymin, xmin, ymax, xmax, clas = boxes[0,i]
    if ymin==-1:
        print("no box")
        continue
    ymin, ymax = lat_[ymin], lat_[ymax]
    xmin, xmax = lon_[xmin], lon_[xmax]
    plt.scatter(np.array([xmin, xmin, xmax, xmax]),np.array([ymin, ymax, ymin, ymax]), c=color[clas])
    rect = patches.Rectangle((xmin,ymin),xmax-xmin,ymax-ymin,linewidth=2,edgecolor=color[clas],facecolor='none')
    ax.add_patch(rect)
from scipy.interpolate import griddata, RectBivariateSpline, RegularGridInterpolator, LinearNDInterpolator, interp2d
from scipy.interpolate import griddata#, RectBivariateSpline, RegularGridInterpolator
Nside = 64
pix = np.arange(12*Nside**2)
coords_hp = hp.pix2ang(Nside, pix, nest=True, lonlat=True)
coords_hp = np.asarray(coords_hp).T
# lon_rad, lat_rad = np.deg2rad(lon), np.deg2rad(lat)
coords_map = hp.ang2vec(lon, lat, lonlat=True).reshape((-1, 3))
coords_map = np.stack([lon, lat], axis=-1).reshape((-1, 2))
# map_hp = griddata(coords_map, images[0,0].flatten(), coords_hp, 'linear')
t = time.time()
map_hp1 = griddata(coords_map, data[0].flatten(), coords_hp, 'linear')
print("time taken:", time.time()-t)

# t = time.time()
# map_hp2 = griddata(coords_map, images[0,0].flatten(), np.asarray(coords_hp).T, 'linear')
# print("time taken:", time.time()-t)

# t = time.time()
# f = interp2d(lon_, lat_, images[0,0])
# map_hp3 = np.diag(f(*coords_hp))  # False result???
# print("time taken:", time.time()-t)

t = time.time()
f = RegularGridInterpolator((lon_, lat_), data[0].T)
map_hp3 = f(coords_hp)
print("time taken:", time.time()-t)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.set_global()
# ax.stock_img()
ax.coastlines()

plt.scatter(coords_hp[:,0], coords_hp[:,1], s=10,
            c=map_hp1, cmap=plt.get_cmap('RdYlBu_r'), alpha=1)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.set_global()
# ax.stock_img()
ax.coastlines()

plt.scatter(coords_hp[:,0], coords_hp[:,1], s=10,
            c=map_hp3, cmap=plt.get_cmap('RdYlBu_r'), alpha=1)
# indic = np.where(~np.isclose(map_hp1, map_hp3, atol=1e-6))
# plt.scatter(coords_hp[indic,0], coords_hp[indic,1], s=10,
#             c='r')
data_path = '../../../data/ExtremeWeather/'
measures, channels, lat_x, lon_x = images.shape
nfeat = 5
Nside = [32, 64]
for nside in Nside:
    npix = hp.nside2npix(nside)
    data = np.empty((measures, npix, channels))
    labels = np.zeros((measures, npix))
    pix = np.arange(npix)
    coords_hp = hp.pix2ang(nside, pix, nest=True, lonlat=True)
    coords_hp = np.asarray(coords_hp).T
    for measure in tqdm(range(measures)):
        for channel in tqdm(range(channels)):
#             data[measure,:,channel] = griddata(coords_map, images[measure,channel].flatten(), coords_hp, 'linear')
            f = RegularGridInterpolator((lon_, lat_), images[measure,channel].T)
            data[measure,:,channel] = f(coords_hp)
        for box in range(boxes.shape[1]):
            ymin, xmin, ymax, xmax, clas = boxes[measure,box]
            if ymin==-1:
#                 print("no box")
                continue
            ymin, ymax = lat_[ymin%lat_x], lat_[ymax%lat_x]
            xmin, xmax = lon_[xmin%lon_x], lon_[xmax%lon_x]
            if xmax>xmin and ymax>ymin:
                indexes = np.where(np.logical_and(np.logical_and(coords_hp[:,0]>=xmin, coords_hp[:,0]<=xmax), 
                                        np.logical_and(coords_hp[:,1]>=ymin, coords_hp[:,1]<=ymax)))
            else:
                indexes = np.where(np.logical_and(np.logical_or(coords_hp[:,0]>=xmin, coords_hp[:,0]<=xmax), 
                                        np.logical_and(coords_hp[:,1]>=ymin, coords_hp[:,1]<=ymax)))
            labels[measure, indexes] = clas + 1
#             if xmax<=xmin or ymax<=ymin:
#                 print(xmin, ymin, xmax, ymax)
#                 print(labels[measure, indexes])
#                 fig = plt.figure(figsize=(25, 50))
#                 ax = fig.add_subplot(211)
#                 ax.scatter(coords_hp[:,0], coords_hp[:,1], 
#                            c=labels[measure,:], cmap=plt.get_cmap('RdYlBu_r'), alpha=1)
#                 ax2 = fig.add_subplot(212, projection=ccrs.PlateCarree())
#                 ax2.set_global()
#                 plt.scatter(coords_hp[:,0], coords_hp[:,1], s=10,
#                     c=labels[measure,:], cmap=plt.get_cmap('RdYlBu_r'), alpha=1)
#                 break
#         if (xmax<=xmin or ymax<=ymin) and ymin!=-1:
#             break
#     break
    file = data_path + 'EW_{}nside_{}'.format(nside, year)
    np.savez(file, data=data, labels=labels)
file = data_path + 'EW_{}nside_{}'.format(Nside, year)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.set_global()
# ax.stock_img()
ax.coastlines()

plt.scatter(coords_hp[:,0], coords_hp[:,1], s=10,
            c=data[0,:,1], cmap=plt.get_cmap('RdYlBu_r'), alpha=1)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.set_global()
# ax.stock_img()
ax.coastlines()

plt.scatter(coords_hp[:,0], coords_hp[:,1], s=10,
            c=labels[0,:], cmap=plt.get_cmap('RdYlBu_r'), alpha=1)
fig = plt.figure(figsize=(25, 25))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(0, 0))
ax.set_global()
ax.coastlines(linewidth=2)

# zmin, zmax = -20, 40

plt.scatter(coords_hp[:,0], coords_hp[:,1], s=100,
            c=map_hp3, cmap=plt.get_cmap('RdYlBu_r'), alpha=1, transform=ccrs.PlateCarree())
def box_labels(predictions):
    # modify predictions to have bounding box
    return predictions

def iou_score(labels, predictions, nclass):
    iou = []
    for i in range(nclass):
        intersect = ((labels==i) + (predictions==i)).eq(2).sum().item()
        union = ((labels==i) + (predictions==i)).sum()
        iou.append(intersect/union)
    return np.asarray(iou)

# def average_precision():
#     pass

# def accuracy():
#     pass<jupyter_output><empty_output><jupyter_text>data<jupyter_code>years = [1979, 1981, 1984, 2005]
path = '../../../data/ExtremeWeather/EW_{}nside_{}.npz'
Nside = 32
datas = []
labels = []
for year in years:
    datapath = path.format(Nside, year)
    file = np.load(datapath)
    data = file['data']  # N x M x F
    label = file['labels']  # N x M
    datas.append(data)
    labels.append(label)
datas = np.stack(datas)
labels = np.stack(labels)
limit=1200
x_train = data[:limit,:,:]
labels_train = labels[:limit,:]
x_val = data[limit:,:,:]
labels_val = labels[limit:,:]

training = LabeledDataset(x_train, labels_train)
validation = LabeledDataset(x_val, labels_val)<jupyter_output><empty_output><jupyter_text>Training<jupyter_code>EXP_NAME = 'TestExtremeWeather'
import tensorflow as tf
params = {'nsides': [nside, nside,nside//2,nside//4,nside//8,nside//16,nside//32, nside//32],
          'F': [20, 20, 50, 80, 100, 200, 400],#np.max(labels_train).astype(int)+1],
          'K': [4]*7,
          'batch_norm': [True]*7}
params['dir_name'] = EXP_NAME
params['num_feat_in'] = x_train.shape[-1] # 2*days_pred+3
params['conv'] = 'chebyshev5'
params['pool'] = 'average'
params['activation'] = 'relu'
params['statistics'] = None#'mean'
params['regularization'] = 0
params['dropout'] = 1
params['num_epochs'] = 1000  # Number of passes through the training data.
params['batch_size'] = 32
params['scheduler'] = lambda step: tf.train.exponential_decay(1e-3, step, decay_steps=2000, decay_rate=1)
#params['optimizer'] = lambda lr: tf.train.GradientDescentOptimizer(lr)
params['optimizer'] = lambda lr: tf.train.RMSPropOptimizer(lr, decay=0.9, momentum=0.)
n_evaluations = 1000
params['eval_frequency'] = int(params['num_epochs'] * (training.N) / params['batch_size'] / n_evaluations)
params['M'] = []
params['Fseg'] = np.max(labels_train).astype(int)+1
print([12*nside**2 for nside in params['nsides']])
model = models.deepsphere(**params)
model.fit(training, validation)
x_train.shape
predictions = model.predict(x_train)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.set_global()
# ax.stock_img()
ax.coastlines()

plt.scatter(coords_hp[:,0], coords_hp[:,1], s=10,
            c=predictions[0,:], cmap=plt.get_cmap('RdYlBu_r'), alpha=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16385246674841997, 'fraction_numerical': 0.03480803697764362, 'mean_word_length': 4.1192660550458715, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '186597', 'license': 'MIT', 'n_tokens_mistral': 4284, 'n_tokens_neox': 4054, 'n_words': 815, 'path': 'Experiments/ExtremeWeather/ExtremeWeatherTest.ipynb', 'repo_name': 'Droxef/PDMdeepsphere', 'size': '10601'}