UP fix notebook message

cf621e6f · BRAUX Emmanuel · f08d2734 · cf621e6f
Commit cf621e6f authored 7 months ago by BRAUX Emmanuel
--- a/turbidity_output_analysis.ipynb
+++ b/turbidity_output_analysis.ipynb
@@ -47,7 +47,7 @@
   "id": "30a87915",
   "metadata": {},
   "source": [
-    "Process Datas (It could take minutes, be patient ...)"
+    "Process Datas : "
   ]
  },
  {
@@ -61,6 +61,16 @@
    "!python main.py xp=base hydra.run.dir=./outputs/main "
   ]
  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a15d307",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    " ... It could take minutes, be patient ..."
+   ]
+  },
  {
   "cell_type": "code",
   "execution_count": null,

 %% Cell type:markdown id:a6c07fb9-260d-49a3-9156-9d27a3ccf1ac tags:

 We perform the inference of a pretrained 4dvarnet model and perform the comparison
 with the DInEOF and eDInEOF models.

 We then choose a range of days specified by a start_date and end_date over which we perform inference.

 %% Cell type:markdown id:f98bbe84 tags:

 Download Datas :
 -  `data/Obs_SPM_log10_aNam.nc`,
 -  `data/Obs_SPM_log10_aNam_removed_50percent_patch_again.nc`,
 -  `data/land_mask_OSE.nc`
 -  `data/checkpoints/val_mse=24.5515-epoch=197.ckpt`

 %% Cell type:code id:e464928f tags:

 ``` python
 !mkdir -p data
 # GT
 !wget -q --show-progress 'https://s3.eu-west-2.wasabisys.com/4dvarnet/german-wadden-sea/Obs_SPM_log10_aNam.nc' -O 'data/Obs_SPM_log10_aNam.n'c
 # patch
 !wget -q --show-progress 'https://s3.eu-west-2.wasabisys.com/4dvarnet/german-wadden-sea/Obs_SPM_log10_aNam_removed_50percent_patch_again.nc' -O 'data/Obs_SPM_log10_aNam_removed_50percent_patch_again.nc'
 # Land mask
 !wget -q --show-progress 'https://s3.eu-west-2.wasabisys.com/4dvarnet/german-wadden-sea/land_mask_OSE.nc' -O 'data/land_mask_OSE.nc'
 # checkpoint
 !mkdir -p data/checkpoints
 !wget -q --show-progress 'https://s3.eu-west-2.wasabisys.com/4dvarnet/german-wadden-sea/val_mse=24.5515-epoch=197.ckpt' -O 'data/checkpoints/val_mse=24.5515-epoch=197.ckpt'
 ```

 %% Cell type:markdown id:30a87915 tags:

-Process Datas (It could take minutes, be patient ...)
+Process Datas :

 %% Cell type:code id:9430a338 tags:

 ``` python
 !rm -rf outputs/main
 !python main.py xp=base hydra.run.dir=./outputs/main
 ```

+%% Cell type:code id:9a15d307 tags:
+
+``` python
+ ... It could take minutes, be patient ...
+```
+
 %% Cell type:code id:de503e67 tags:

 ``` python
 # We import the libraries needed for the pipeline

 import xarray as xr
 import matplotlib.pyplot as plt
 import numpy as np
 import scipy.sparse.linalg as sp
 import matplotlib.animation as animation
 import warnings
 plt.rcParams["animation.html"] = "jshtml"
 plt.rcParams['figure.dpi'] = 150
 plt.ioff()
 ```

 %% Cell type:code id:4fa7b0dd tags:

 ``` python
 # The two performance metrics we use for cross-validation in (e)DInEOF algorithms)

 def rmse(x,y):
    return np.sqrt(np.mean((x-y)**2))

 def RE(gt,pred):
    ep=0.0001
    return np.mean(np.abs((10**gt+ep)-(10**pred+ep))/np.abs(10**gt+ep))*100
 ```

 %% Cell type:code id:b132951b tags:

 ``` python
 # DInEOF algorithm (from https://doi.org/10.1016/j.ocemod.2004.08.001)

 def DInEOF(X, mask_miss, mask_test=None, k=1, iter_=1000):
    X_rec_temp = np.zeros(X.shape)
    score = np.zeros([3,iter_])
    for i in range(iter_):
        X_rec_out = X_rec_temp.copy()
        X_rec_temp[~mask_miss] = X[~mask_miss].copy()
        #u, s, vh = sp.svds(X_rec_temp,k=k)
        #X_rec_temp = u@np.diag(s)@vh
        A = X_rec_temp.T@X_rec_temp
        p2, v = sp.eigs(A,k)
        if np.sum(np.imag(p2))!=0:
            warning.warn("Caution non-zero imaginary part has been discard")
        p = np.sqrt(np.real(p2))
        v = np.real(v)
        u = X_rec_temp@(v@np.diag(1./p))
        X_rec_temp = u@np.diag(p)@v.T
        score[0,i] = i
        score[1,i] = k
        score[2,i] = rmse(X[mask_test], X_rec_temp[mask_test])
        if i>0:
            if score[2,i]>score[2,i-1]:
                score = score[:,0:i+1]
                return X_rec_out, score
            elif np.abs(score[2,i]-score[2,i-1])<1.e-4:
                k+=1
    X_rec_out = X_rec_temp
    return X_rec_out, score
 ```

 %% Cell type:code id:2857040f tags:

 ``` python
 # eDInEOF algorithm (from https://doi.org/10.5194/os-5-475-2009)

 def eDInEOF(X, mask_miss, mask_test, dt=1, alpha=1.e-2 , s_=5, k=1, iter_=1000):
    N = X.shape[1]
    F = -2*np.eye(N,N)+np.eye(N,N,-1)+np.eye(N,N,1)
    F[1,0] = 2
    F[-2,-1] = 2
    F = (alpha/dt**2)*F

    X_rec_temp = np.zeros(X.shape)
    score = np.zeros([3,iter_])
    for i in range(iter_):
        X_rec_out = X_rec_temp.copy()
        X_rec_temp[~mask_miss] = X[~mask_miss].copy()
        for s in range(s_):
            X_rec_temp = X_rec_temp+X_rec_temp@F
        #u, sigma, vh = sp.svds(X_rec_temp, k=k)
        #X_rec_temp = u@np.diag(sigma)@vh
        A = X_rec_temp.T@X_rec_temp
        p2, v = sp.eigs(A,k)
        if np.sum(np.imag(p2))!=0:
            warning.warn("Caution non-zero imaginary part has been discard")
        p = np.sqrt(np.real(p2))
        v = np.real(v)
        u = X_rec_temp@(v@np.diag(1./p))
        X_rec_temp = u@np.diag(p)@v.T
        score[0,i] = i
        score[1,i] = k
        score[2,i] = rmse(X[mask_test], X_rec_temp[mask_test])
        if i>0:
            if score[2,i]>score[2,i-1]:
                score = score[:,0:i+1]
                return X_rec_out, score
            elif np.abs(score[2,i]-score[2,i-1])<1.e-4:
                k+=1
    X_rec_out = X_rec_temp
    return X_rec_out, score

 ```

 %% Cell type:code id:4cb86eae tags:

 ``` python
 #This cell is to prepare for computing DInEOF and eDInEOF
 var = 'SPM'

 start_date='2020-06-01'
 end_date='2020-06-26'

 #OSE
 data = xr.open_dataset('data/Obs_SPM_log10_aNam_removed_50percent_patch_again.nc').sel(time=slice(start_date,end_date)) #Gappy observations
 land_mask = xr.open_dataset('data/land_mask_OSE.nc')


 nlat = data.lat.size
 nlon = data.lon.size

 nb_sea_pix = np.sum(np.isnan(land_mask.mask.values)) # Nb of sea pixels to be filled
 idx_sea = np.isnan(land_mask.mask.values.flatten()) # Sea pixels locations

 dt = len(data.time) # Length of the time-serie
 X = np.full([nb_sea_pix,dt],np.nan) # Init. of the matrix to store vectorized sea-pixels

 # Loop to build the data matrix to be used in completion algorithms
 for t in range(dt):
    data_t = data[var][t].values.flatten()
    X[:,t] = data_t[idx_sea]

 mask_mv = np.isnan(X) # Missing data mask matrix
 mask_data = 1-mask_mv # Data mask matrix
 data_index = np.where(mask_data) # Data index

 mask_crossval = np.full([nb_sea_pix,dt], False, dtype=bool) # Init. mask matrix for cross-validation data
 rand_pick = np.random.permutation(data_index[0].size)[0:int(.01*data_index[0].size)]
 data_pick = (data_index[0][rand_pick],data_index[1][rand_pick]) # Random selection of data for cros-validation (1%)
 mask_crossval[data_pick] = True

 mask_mv[mask_crossval] = True

 # Data normalization
 mu = np.nanmean(X)
 std = np.nanstd(X)**2
 X = (X-mu)/std
 X = np.where(np.isnan(X),0,X)
 ```

 %% Cell type:code id:ae7da4dd tags:

 ``` python
 # Apply DInEOF to the selected data
 X_rec, score = DInEOF(X, mask_mv, mask_crossval)
 X_rec = X_rec*std+mu

 # Reshaping
 rec_DInEOF = np.zeros([dt,nlat,nlon])
 for t in range(dt):
    data_t = np.full([nlat*nlon], np.nan)
    data_t[idx_sea] = X_rec[:,t]
    rec_DInEOF[t] = data_t.reshape(nlat,nlon)

 #saving reconstruction
 rec_DInEOF = xr.DataArray(rec_DInEOF,name=var,coords=[data.time[0:dt],data.lat,data.lon],dims=['time','lat','lon'])

 #OSE
 rec_DInEOF.to_netcdf('data/DInEOF_log10.nc')
 ```

 %% Cell type:code id:f52da7d4-0622-41ec-a00b-51551518fa95 tags:

 ``` python
 # Apply eDInEOF to the selected data
 X_rec, score_e = eDInEOF(X, mask_mv, mask_crossval)
 X_rec = X_rec*std+mu

 # Reshaping
 rec_eDInEOF = np.zeros([dt,nlat,nlon])
 for t in range(dt):
    data_t = np.full([nlat*nlon], np.nan)
    data_t[idx_sea] = X_rec[:,t]
    rec_eDInEOF[t] = data_t.reshape(nlat,nlon)

 # saving reconstruction
 rec_eDInEOF = xr.DataArray(rec_eDInEOF,name=var,coords=[data.time[0:dt],data.lat,data.lon],dims=['time','lat','lon'])

 #OSE
 rec_eDInEOF.to_netcdf('data/eDInEOF_log10.nc')
 ```

 %% Cell type:code id:63c77a1b tags:

 ``` python
 # Compute the RMSE, RE of the three algorithms

 # #load data OSE
 rec_DInEOF = xr.open_dataset('data/DInEOF_log10.nc')#.sel(time="2018")
 rec_eDInEOF = xr.open_dataset('data/eDInEOF_log10.nc')#.sel(time="2018")
 data4DVarnet = xr.open_dataset("outputs/main/base/TestonOSE_DutchWaddenSea/test_data.nc")

 GT = xr.open_dataset('data/Obs_SPM_log10_aNam.nc').sel(time=slice(start_date,end_date))
 Obs = xr.open_dataset('data/Obs_SPM_log10_aNam_removed_50percent_patch_again.nc').sel(time=slice(start_date,end_date))


 #declare the pixels in which we process RMSE and RE
 mask_obs_minus_GT = ~np.isnan(GT.SPM.values) & np.isnan(Obs.SPM)#mask_obs_minus_GT = True if Obs is nan and GT is not nan (and False otherwise)

 #score_DInEOF
 score_DInEOF = rmse(GT.SPM.values[mask_obs_minus_GT],rec_DInEOF.SPM.values[mask_obs_minus_GT])
 print("RMSE DinEOF: ",score_DInEOF)

 #RE_DInEOF
 RE_DInEOF = RE(GT.SPM.values[mask_obs_minus_GT],rec_DInEOF.SPM.values[mask_obs_minus_GT])
 print("RE DinEOF: ", RE_DInEOF)

 #score_eDInEOF
 score_eDInEOF = rmse(GT.SPM.values[mask_obs_minus_GT],rec_eDInEOF.SPM.values[mask_obs_minus_GT])
 print("RMSE eDinEOF: ",score_eDInEOF)

 #RE_eDInEOF
 RE_eDInEOF = RE(GT.SPM.values[mask_obs_minus_GT],rec_eDInEOF.SPM.values[mask_obs_minus_GT])
 print("RE eDinEOF: ",RE_eDInEOF)

 #score_4DVarNetwith
 score_4DVarNet = rmse(GT.SPM.values[mask_obs_minus_GT],data4DVarnet.out.values[mask_obs_minus_GT])
 print("RMSE 4DVarNEt: ",score_4DVarNet)

 #RE_4DVarNet
 RE_4DVarNet = RE(data4DVarnet.tgt.values[mask_obs_minus_GT],data4DVarnet.out.values[mask_obs_minus_GT])
 print("RE 4DVarNEt: ",RE_4DVarNet)
 ```

 %% Cell type:code id:402d80d2 tags:

 ``` python
 #We load the land mask and the satelite observations

 land_mask=xr.open_dataset('data/land_mask_OSE.nc')
 GT = xr.open_dataset('data/Obs_SPM_log10_aNam.nc').sel(time=slice(start_date,end_date))
 ```

 %% Cell type:code id:e982290b-5322-49f3-a586-db3c47034893 tags:

 ``` python
 #plot all algorithms for some specific days
 rang_min=5#220,130
 rang_max=rang_min+1
 step=1

 lon = GT.lon
 lat = GT.lat

 nlat=lat.size
 nlon=lon.size

 # transpose land_mask
 land_mask_transposed = land_mask.transpose('lat', 'lon')


 v_min=0.0#np.min(data.tgt[rang_min:rang_max:10])
 v_max=2.0#np.max(data.tgt[rang_min:rang_max:10])

 plots_path='plots'
 import os
 if not os.path.exists(plots_path):
    os.makedirs(plots_path)

 for t in range(rang_min, rang_max, step):
    fig, axes = plt.subplots(2, 3, figsize=(4 * 4, 2 * 4 * (nlat / nlon)))  # Adjust size as needed for 2 rows, 3 columns


    # First subplot (Row 1, Column 1)
    mappable0 = axes[0, 0].pcolormesh(lon, lat, data4DVarnet.tgt[t], cmap='jet', vmin=v_min, vmax=v_max)
    axes[0, 0].pcolormesh(lon, lat, land_mask.mask, cmap='Greys', vmin=0, vmax=2)
    axes[0, 0].set_title(GT.time[t].dt.strftime("%B %d, %Y").values + " - Gappy Target")


    # Second subplot (Row 1, Column 2)
    mappable1 = axes[0, 1].pcolormesh(lon, lat, data4DVarnet.inp[t], cmap='jet', vmin=v_min, vmax=v_max)
    axes[0, 1].pcolormesh(lon, lat, land_mask.mask, cmap='Greys', vmin=0, vmax=2)
    axes[0, 1].set_title("Simulated Observation")


    # Third subplot (Row 1, Column 3)
    # If you have a third type of data for this row, use it here
    # Example:
    mappable2 = axes[0, 2].pcolormesh(lon, lat, rec_DInEOF.SPM.values[t], cmap='jet', vmin=v_min, vmax=v_max)
    axes[0, 1].pcolormesh(lon, lat, land_mask.mask, cmap='Greys', vmin=0, vmax=2)
    axes[0, 2].set_title("DInEOF")


    # Fourth subplot (Row 2, Column 1)
    mappable3 = axes[1, 0].pcolormesh(lon, lat, rec_eDInEOF.SPM.values[t], cmap='jet', vmin=v_min, vmax=v_max)
    axes[1, 0].pcolormesh(lon, lat, land_mask.mask, cmap='Greys', vmin=0, vmax=2)
    axes[1, 0].set_title("eDInEOF")


    # Fifth subplot (Row 2, Column 2)
    mappable4 = axes[1, 1].pcolormesh(lon, lat, data4DVarnet.out[t], cmap='jet', vmin=v_min, vmax=v_max)
    axes[1, 1].pcolormesh(lon, lat, land_mask.mask, cmap='Greys', vmin=0, vmax=2)
    axes[1, 1].set_title("4DVarNEt")


   # Sixth subplot (Row 2, Column 3)
   #mappable5 = axes[1, 2].pcolormesh(lon, lat, data4DVarnet_Conv3D_OSSE2OSE.out[t], cmap='jet', vmin=v_min, vmax=v_max)
   #axes[1, 2].pcolormesh(lon, lat, land_mask.mask, cmap='Greys', vmin=0, vmax=2)
   #axes[1, 2].set_title("4DVarNEt_Conv3D_OSSE2OSE")


    # Colorbar setup
    fig.subplots_adjust(right=0.8)
    cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.7])  # Position for the colorbar
    fig.colorbar(mappable0, cax=cbar_ax, extend='min')


    # Save the figure
    fig.savefig(f"plots/plot_4DVarNet_vs_Unet_vs_DinEOF_SCHISM_time{t}.png", dpi=300)


    plt.show()
 ```