Case study of prediction performance in Pearl River Estuary

[1]:

import h5py
import scipy
import numpy as np
import torch

[2]:

def load(path):
    path = path + "/with_imputation"
    preds = []
    for i in range(10):
        preds.append(np.load(path+"/prediction_{}.npy".format(i)))
    return np.stack(preds, axis=1)

[3]:

base_dir = "../log/prediction/PRE/"
prediction_xg_wo = torch.from_numpy(np.load(base_dir+"XGBoost/without_imputation/prediction_0.npy", allow_pickle=True))
prediction_our = torch.from_numpy(load(base_dir+"STIMP"))
prediction_predrnn_wo = torch.from_numpy(np.load(base_dir+"PredRNN/without_imputation/prediction_0.npy", allow_pickle=True))
prediction_pde = np.load("../data/PRE/cmoms.npy")

[4]:

label = np.load("../data/PRE/trues.npy")
label_masks = np.load("../data/PRE/true_masks.npy")
index = [46*i for i in range(306//46)]
label_masks = label_masks.squeeze()
label = label.squeeze()
label = torch.from_numpy(label)
label_masks = torch.from_numpy(label_masks)

Fig 3i: Predicted Chl_a of CMOMS, STIMP and PredRNN at Position 1

[6]:

#CMOMS
from copy import deepcopy
import seaborn as sns
import matplotlib.pyplot as plt
n=803

index = [46*i for i in range(306//46)]
tmp  = deepcopy(label[index].reshape(276,-1))
tmp_mask  = deepcopy(label_masks[index].reshape(276,-1))
tmp[~tmp_mask.bool()]=np.nan
plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_pde[:276])
predict = predict[:,n]

plt.plot(np.arange(276), predict, label="CMOMS")

plt.ylim(-0.5,1.5)
plt.xticks([])
plt.yticks([-0.4,0.5,1.4], fontsize=16)
plt.ylabel("Chl_a (Log(ug/L))", fontsize=20)
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_6_0.png

[7]:

#STIMP
plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")
predict = deepcopy(prediction_our[index].transpose(1,2).reshape(276,10,4443))
predict = predict[:,:,n]
mean = predict.mean(1)
std = predict.std(1)

plt.plot(np.arange(276), mean, label="STImp")
plt.fill_between(np.arange(276), mean-std, mean+std, alpha=0.3)

plt.xticks(fontsize=14)
plt.yticks([])
plt.ylim(-0.5,1.5)
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_7_0.png

[8]:

#PredRNN
plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")
predict = deepcopy(prediction_predrnn_wo[index].reshape(276,4443))
predict = predict[:,n]
plt.plot(np.arange(276), predict, label="STImp")

plt.ylim(-0.5,1.5)
plt.yticks([])
plt.xticks([])
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_8_0.png

Fig 3i: Predicted Chl_a of CMOMS, STIMP and PredRNN at Position 2

[9]:

n=4006

index = [46*i for i in range(306//46)]
tmp  = deepcopy(label[index].reshape(276,-1))
tmp_mask  = deepcopy(label_masks[index].reshape(276,-1))
tmp[~tmp_mask.bool()]=np.nan
plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_pde[:276])
predict = predict[:,n]

plt.plot(np.arange(276), predict, label="CMOMS")

plt.ylim(-2.8,0.1)
plt.xticks([])
plt.yticks(np.arange(-2.5, 0.5, 1), fontsize=16)
plt.ylabel("Chl_a (Log(ug/L))", fontsize=20)
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_10_0.png

[10]:

#STIMP
plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_our[index].transpose(1,2).reshape(276,10,4443))
predict = predict[:,:,n]
mean = predict.mean(1)
std = predict.std(1)

plt.plot(np.arange(276), mean, label="STImp")
plt.fill_between(np.arange(276), mean-std, mean+std, alpha=0.3)

# plt.legend()
plt.xticks([])
plt.yticks([])
plt.ylim(-2.8,0.1)
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_11_0.png

[11]:

plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_predrnn_wo[index].reshape(276,4443))
predict = predict[:,n]

plt.plot(np.arange(276), predict, label="STImp")

# plt.legend()
plt.ylim(-2.8,0.1)
plt.yticks([])
plt.xticks([])
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_12_0.png

Fig 3i: Predicted Chl_a of CMOMS, STIMP and PredRNN at Position 3

[12]:

n=2667

index = [46*i for i in range(306//46)]
tmp  = deepcopy(label[index].reshape(276,-1))
tmp_mask  = deepcopy(label_masks[index].reshape(276,-1))
tmp[~tmp_mask.bool()]=np.nan
plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_pde[:276])
predict = predict[:,n]

plt.plot(np.arange(276), predict, label="CMOMS")

plt.ylim(-2.2,0.5)
plt.xticks([])
plt.yticks([-2.1,-1.1,0.4], fontsize=16)
plt.ylabel("Chl_a (Log(ug/L))", fontsize=20)
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_14_0.png

[13]:

plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_our[index].transpose(1,2).reshape(276,10,4443))
predict = predict[:,:,n]
mean = predict.mean(1)
std = predict.std(1)

plt.plot(np.arange(276), mean, label="STImp")
plt.fill_between(np.arange(276), mean-std, mean+std, alpha=0.3)

# plt.legend()
plt.xticks([])
plt.yticks([])
plt.ylim(-2.2,0.5)
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_15_0.png

[14]:

plt.scatter(np.arange(276), tmp[:,n], c='red', s=10,label="Observed")

predict = deepcopy(prediction_predrnn_wo[index].reshape(276,4443))
predict = predict[:,n]

plt.plot(np.arange(276), predict, label="STImp")

# plt.legend()
plt.ylim(-2.2,0.5)
plt.yticks([])
plt.xticks([])
plt.show()

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_16_0.png

Fig 3j: True Chl_a, predicted Chl_a of STIMP and DINEOF on 2 February 2016

[15]:

from mpl_toolkits import basemap
from sklearn.cluster import KMeans, DBSCAN, SpectralClustering
import cartopy.crs as ccrs
from copy import deepcopy
import h5py
from numpy import meshgrid
import numpy as np

t = 43
is_sea = np.load("../data/PRE/is_sea.npy")
tmp = deepcopy(label[46,t,:])
tmp[~label_masks[0,t,:].bool()] = np.nan
ob = np.zeros((60,96))
ob[is_sea.astype(bool)] = tmp
ob[~is_sea.astype(bool)]= np.nan
lon = np.load("../data/PRE/lon.npy")
lati = np.load("../data/PRE/lati.npy")

lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, ob, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend="both")

[15]:

<matplotlib.contour.QuadContourSet at 0x75e3c028ac40>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_18_1.png

Fig 3j: Predicted Chl_a of STIMP on 2 February 2016

[16]:

tmp = deepcopy(prediction_our[46,:, t,:].mean(0))
pre = np.zeros((60,96))
pre[is_sea.astype(bool)] = tmp
pre[~is_sea.astype(bool)]= np.nan


lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, pre, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend='both')

[16]:

<matplotlib.contour.QuadContourSet at 0x75e6ab736670>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_20_1.png

[18]:

tmp = deepcopy(prediction_pde[46+t,:])
pre = np.zeros((60,96))
pre[is_sea.astype(bool)] = tmp
pre[~is_sea.astype(bool)]= np.nan


lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, pre, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend='both')

# map.contourf(x, y, tmp2, levels=np.linspace(-1.5, 1.5, 40),cmap="Greys")
map.colorbar(boundaries=np.linspace(-1.5, 1.5, 20), ticks=np.linspace(-1.5, 1.5, 6))

[18]:

<matplotlib.colorbar.Colorbar at 0x75e3ba9111f0>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_21_1.png

Fig 3j: True Chl_a, predicted Chl_a of STIMP and DINEOF on 16 July 2016

[19]:

t = 0
is_sea = np.load("../data/PRE/is_sea.npy")
tmp = deepcopy(label[46,t,:])
tmp[~label_masks[0,t,:].bool()] = np.nan
ob = np.zeros((60,96))
ob[is_sea.astype(bool)] = tmp
ob[~is_sea.astype(bool)]= np.nan
lon = np.load("../data/PRE/lon.npy")
lati = np.load("../data/PRE/lati.npy")

lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, ob, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend="both")

[19]:

<matplotlib.contour.QuadContourSet at 0x75e3bb3ec700>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_23_1.png

[20]:

tmp = deepcopy(prediction_our[46,:, t,:].mean(0))
pre = np.zeros((60,96))
pre[is_sea.astype(bool)] = tmp
pre[~is_sea.astype(bool)]= np.nan


lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, pre, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend='both')

[20]:

<matplotlib.contour.QuadContourSet at 0x75e3bb32d700>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_24_1.png

[21]:

tmp = deepcopy(prediction_pde[46+t,:])
pre = np.zeros((60,96))
pre[is_sea.astype(bool)] = tmp
pre[~is_sea.astype(bool)]= np.nan


lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, pre, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend='both')

# map.contourf(x, y, tmp2, levels=np.linspace(-1.5, 1.5, 40),cmap="Greys")
map.colorbar(boundaries=np.linspace(-1.5, 1.5, 20), ticks=np.linspace(-1.5, 1.5, 6))

[21]:

<matplotlib.colorbar.Colorbar at 0x75e3bb31ab50>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_25_1.png

Fig 3j: True Chl_a, predicted Chl_a of STIMP and DINEOF on 8 January 2016

[23]:

t = 20
is_sea = np.load("../data/PRE/is_sea.npy")
tmp = deepcopy(label[46,t,:])
tmp[~label_masks[0,t,:].bool()] = np.nan
ob = np.zeros((60,96))
ob[is_sea.astype(bool)] = tmp
ob[~is_sea.astype(bool)]= np.nan
lon = np.load("../data/PRE/lon.npy")
lati = np.load("../data/PRE/lati.npy")

lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, ob, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend="both")

[23]:

<matplotlib.contour.QuadContourSet at 0x75e3bb636160>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_27_1.png

[24]:

tmp = deepcopy(prediction_our[46,:, t,:].mean(0))
pre = np.zeros((60,96))
pre[is_sea.astype(bool)] = tmp
pre[~is_sea.astype(bool)]= np.nan


lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, pre, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend='both')

[24]:

<matplotlib.contour.QuadContourSet at 0x75e3bb3d45b0>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_28_1.png

[25]:

tmp = deepcopy(prediction_pde[46+t,:])
pre = np.zeros((60,96))
pre[is_sea.astype(bool)] = tmp
pre[~is_sea.astype(bool)]= np.nan


lon1, lon2, lati1, lati2 = lon.min(), lon.max(), lati.min(), lati.max()
map = basemap.Basemap(llcrnrlon=lon1, llcrnrlat=lati1,urcrnrlon=lon2, urcrnrlat=lati2, projection='cyl', resolution='f')
# map.fillcontinents(color='white')
map.drawlsmask(land_color='white', ocean_color='lightgray', resolution='f',grid=1.25)
# map.bluemarble()
map.drawcoastlines()
map.contourf(lon, lati, pre, levels=np.linspace(-1.5, 1.5, 40),cmap="rainbow",extend='both')

# map.contourf(x, y, tmp2, levels=np.linspace(-1.5, 1.5, 40),cmap="Greys")
map.colorbar(boundaries=np.linspace(-1.5, 1.5, 20), ticks=np.linspace(-1.5, 1.5, 6))

[25]:

<matplotlib.colorbar.Colorbar at 0x75e3bb3cbd60>

../../_images/analysis_PRE_04-prediction-pearl-river-estuary-case-study_29_1.png