Skip to content

Commit

Permalink
Issue #46 use module functions from notebooks
Browse files Browse the repository at this point in the history
  • Loading branch information
@soxofaan
soxofaan committed Apr 25, 2022
1 parent 72f656d commit 5448478
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 225 deletions.
158 changes: 3 additions & 155 deletions src/fusets/UVScripts/AI4FOOD_OpenEO_MOGPR_S1S2.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,8 @@
"from datetime import datetime\n",
"import itertools\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.dates as mdates"
"import matplotlib.dates as mdates\n",
"from fusets.mogrp import MOGRP_GPY_retrieval"
]
},
{
Expand Down Expand Up @@ -302,159 +303,6 @@
"#fig.savefig(\"./fig_out/\"+crop+\"TS.png\", bbox_inches='tight')"
]
},
{
"cell_type": "markdown",
"id": "1ab14c1d",
"metadata": {},
"source": [
"### MOGPR (GPY)\n",
"Function definition"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "8097b5ac",
"metadata": {},
"outputs": [],
"source": [
"def create_empty_2D_list(nrows,ncols):\n",
" \"\"\"\n",
" Function creating an empty list of lists (matrix-like structure) \n",
" \n",
" Args:\n",
" nrows [int] : # of rows of the output list of lists\n",
" ncols [int] : # of cols of the output list of lists\n",
"\n",
" Returns:\n",
" empty_list [list] : empty list of lists (matrix-like structure) \n",
" \n",
" \"\"\"\n",
" empty_list = []\n",
" for j in range(nrows):\n",
" column = []\n",
" for i in range(ncols):\n",
" column.append(0)\n",
" empty_list.append(column)\n",
" return empty_list\n",
"\n",
"\n",
"def MOGRP_GPY_retrieval(data_in, time_in, master_ind, output_timevec, nt):\n",
" \"\"\"\n",
" Function performing the multioutput gaussian-process regression at pixel level for gapfilling purposes\n",
"\n",
" Args:\n",
" data_in [array] : 3D (2DSpace, Time) array containing data to be processed \n",
" time_in [array] : vector containing the dates of each layer in the time dimension\n",
" master_ind [int] : Index identifying the Master output\n",
" output_timevec [array] :vector containing the dates on which output must be estimated\n",
" nt [int] : # of time the GP training must be performed (def=1) \n",
" Returns:\n",
" Out_mean [array] : 3D (2DSpace, Time) mean value of the prediction at pixel level\n",
" Out_std [array] : 3D (2DSpace, Time) standard deviation of the prediction at pixel level\n",
" Out_QFlag [array]: 2D map of Quality Flag for any numerical error in the model determination\n",
" Out_model [list] : Matrix-like structure containing the model information at pixel level \n",
" \"\"\" \n",
" noutput_timeseries = len(data_in)\n",
"\n",
" x_size = data_in[0].shape[1]\n",
" y_size = data_in[0].shape[2]\n",
" imout_sz = (output_timevec.shape[0],x_size,y_size)\n",
"\n",
" Xtest = output_timevec.reshape(output_timevec.shape[0],1)\n",
" out_qflag = np.ones((x_size,y_size), dtype=bool)\n",
" out_mean = []\n",
" out_std = []\n",
" out_model = create_empty_2D_list(x_size,y_size)\n",
"\n",
" for _ in range(noutput_timeseries):\n",
" out_mean.append(np.full(imout_sz,np.nan))\n",
" out_std.append(np.full(imout_sz,np.nan))\n",
" #Variable is initialized to take into account possibility that within one block no valid pixels is present, due to the mask\n",
" for x, y in itertools.product(range(x_size), range(y_size)):\n",
"\n",
" X_vec = []\n",
" Y_vec = []\n",
" Y_mean_vec = []\n",
" Y_std_vec = []\n",
" \n",
" for ind in range(noutput_timeseries):\n",
"\n",
" X_tmp = time_in[ind]\n",
" Y_tmp = data_in[ind][:,x,y] \n",
" X_tmp = X_tmp[~np.isnan(Y_tmp),np.newaxis]\n",
" Y_tmp = Y_tmp[~np.isnan(Y_tmp),np.newaxis]\n",
" X_vec.append(X_tmp)\n",
" Y_vec.append(Y_tmp)\n",
" del X_tmp,Y_tmp\n",
" \n",
" \n",
" if np.size(Y_vec[master_ind]) >0:\n",
" \n",
" # Data Normalization\n",
" for ind in range(noutput_timeseries):\n",
" Y_mean_vec.append(np.mean(Y_vec[ind])) \n",
" Y_std_vec.append(np.std(Y_vec[ind]))\n",
" Y_vec[ind] = (Y_vec[ind]-Y_mean_vec[ind])/Y_std_vec[ind]\n",
" \n",
" # Multi-output train and test sets\n",
" Xtrain = X_vec\n",
" Ytrain = Y_vec\n",
" \n",
" nsamples, npixels = Xtest.shape\n",
" noutputs = len(Ytrain)\n",
" \n",
" for i_test in range(nt):\n",
" Yp = np.zeros((nsamples, noutputs))\n",
" Vp = np.zeros((nsamples, noutputs))\n",
" #K = GPy.kern.RBF(1, lengthscale=32.9169)\n",
" K = GPy.kern.Matern32(1, lengthscale=32.9169) # Use RBF or Matern32 as kernels\n",
" LCM = GPy.util.multioutput.LCM(input_dim = 1, \n",
" num_outputs = noutputs,\n",
" kernels_list = [K] * noutputs, W_rank=1) \n",
" model = GPy.models.GPCoregionalizedRegression(Xtrain, Ytrain, kernel=LCM.copy()) \n",
" model['.*Mat32.var'].constrain_fixed(1.) \n",
" #model['.*rbf.var'].constrain_fixed(1.)\n",
" if not np.isnan(Ytrain[1]).all():\n",
"\n",
" try:\n",
" #if trained_model is None:\n",
" model.optimize()\n",
" list_tmp = [model.param_array]\n",
" \n",
" for _ in range(noutput_timeseries):\n",
" list_tmp.append(eval('model.sum.ICM'+str(_)+'.B.B'))\n",
" out_model[x][y]=list_tmp\n",
" \n",
" \n",
" except:\n",
" out_qflag[x,y]=False\n",
" continue\n",
" \n",
"\n",
" for out in range(noutputs):\n",
" newX = Xtest.copy() \n",
" \n",
" newX = np.hstack([newX, out * np.ones((newX.shape[0], 1))])\n",
" noise_dict = {'output_index': newX[:, -1:].astype(int)}\n",
" Yp[:,None, out],Vp[:,None,out] =model.predict(newX, Y_metadata=noise_dict)\n",
"\n",
" if i_test==0:\n",
" \n",
" for ind in range(noutput_timeseries):\n",
" out_mean[ind][:,None,x,y] = (Yp[:,None, 0]*Y_std_vec[ind]+Y_mean_vec[ind])/Nt \n",
" out_std[ind][:,None,x,y] = (Vp[:,None, 0]*Y_std_vec[ind])/Nt \n",
"\n",
" else:\n",
" for ind in range(noutput_timeseries):\n",
" out_mean[ind][:,None,x,y] = out_mean[ind][:,None,x,y] + (Yp[:,None, 0]*Y_std_vec[ind]+Y_mean_vec[ind])/Nt \n",
" out_std[ind][:,None,x,y] = out_std[ind][:,None,x,y] + (Vp[:,None, 0]*Y_std_vec[ind])/Nt \n",
" \n",
" del Yp,Vp\n",
" \n",
" return out_mean, out_std, out_qflag, out_model\n"
]
},
{
"cell_type": "markdown",
"id": "9caf78c1",
Expand Down Expand Up @@ -942,4 +790,4 @@
},
"nbformat": 4,
"nbformat_minor": 5
}
}
73 changes: 3 additions & 70 deletions src/fusets/UVScripts/AI4FOOD_Whittaker_OpenEO_NDVI_S2.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,8 @@
"import math\n",
"from shapely.geometry import box\n",
"from openeo.rest.conversions import timeseries_json_to_pandas\n",
"import time"
"import time\n",
"from fusets.whittaker import whittaker_f"
]
},
{
Expand All @@ -52,74 +53,6 @@
"connection = openeo.connect(\"openeo.vito.be\").authenticate_oidc()"
]
},
{
"cell_type": "markdown",
"id": "4fffc877",
"metadata": {},
"source": [
"Custom functions"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "37943ddd",
"metadata": {},
"outputs": [],
"source": [
"def whittaker_f(x,y,lmbd,d):\n",
"\n",
" #minimum and maximum dates\n",
" min=x[0].toordinal()\n",
" max=x[-1].toordinal()\n",
" l=max-min\n",
"\n",
" v=np.full(l+1,-3000)\n",
" t=np.full(l+1,0,dtype='float')\n",
"\n",
" D=[]\n",
" for i in x:\n",
" D.append(i.toordinal())\n",
"\n",
"\n",
" D1=np.array(D)-D[0]\n",
" D11=D1[~np.isnan(y)]\n",
" v[D11]=1\n",
" t[D11]=y[~np.isnan(y)]\n",
"\n",
" # dates\n",
" xx = [x[0] + timedelta(days=i) for i in range(l+1)] #try using pandas instead?\n",
"\n",
" # create weights\n",
" w = np.array((v!=-3000)*1,dtype='double')\n",
"\n",
" # apply filter\n",
" z_ = ws2d(t,lmbd,w)\n",
" z1_ = np.array(z_)\n",
" \n",
" #return z1_,xx\n",
" \n",
" if isinstance(d, int): \n",
" if d>0 and d<z1_.size: \n",
" n=z1_.size\n",
" ind=[i*d for i in range(math.ceil(n/d))]\n",
" Zd=z1_[ind]\n",
" XXd=[xx[ii] for ii in ind] \n",
" else:\n",
" Zd=[]\n",
" XXd=[]\n",
" print (\"d must be positive and smaller than the total number of dates\")\n",
" else:\n",
" Zd=[]\n",
" XXd=[]\n",
" print (\"d must be an integer\")\n",
" \n",
" return z1_,xx,Zd,XXd\n",
" \n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 4,
Expand Down Expand Up @@ -399,4 +332,4 @@
},
"nbformat": 4,
"nbformat_minor": 5
}
}

0 comments on commit 5448478

Please sign in to comment.