lr.py 2.30 KiB
import pandas as pd 
import numpy as np 
from sklearn.preprocessing import StandardScaler
from sklearn import linear_model
from sklearn import metrics
import seaborn as sns
import matplotlib.pyplot as plt

def simple_lr (X,y,n_timestamps,bands):
    n_bands = len(bands) #int(X.shape[1]/n_timestamps)
    times = pd.date_range(start="20170502", end = "20171030", freq='5D')
    for i in range(n_bands):
        dic = {}
        for j in range (n_timestamps):
            x = X[:,n_bands*j+i]
            x = x.reshape(-1,1)
            lm = linear_model.LinearRegression()
            lm.fit(x,y)

            y_predict = lm.predict(x)

            dic.setdefault("estimator",[]).append("%s_%s"%(bands[i],times.strftime('%m-%d')[j]))
            dic.setdefault("r2",[]).append(metrics.r2_score(y,y_predict))
            dic.setdefault("mae",[]).append(metrics.mean_absolute_error(y,y_predict))
            dic.setdefault("mse",[]).append(metrics.mean_squared_error(y,y_predict))
            dic.setdefault("rmse",[]).append(np.sqrt(metrics.mean_squared_error(y,y_predict)))
    
        df = pd.DataFrame.from_dict(dic)
        print (df.loc[df[['r2']].idxmax()])
        # print(df.loc[df[['mae','mse','rmse']].idxmin()])


if __name__ == '__main__' :
    radar_X1 = np.load("./data/niakhar_2017_rad_seq.npy")
    opt_X1 =  np.load("./data/niakhar_2017_opt_seq.npy")
    indices_X1 =  np.load("./data/niakhar_2017_indices_seq.npy")
    y1 =  np.load("./data/niakhar_2017_yields.npy")

    radar_X2 =  np.load("./data/niakhar-simco_2018_rad_seq.npy")
    opt_X2 =  np.load("./data/niakhar-simco_2018_opt_seq.npy")
    indices_X2 =  np.load("./data/niakhar-simco_2018_indices_seq.npy")
    y2 =  np.load("./data/niakhar-simco_2018_yields.npy")

    radar_X = np.vstack((radar_X1,radar_X2))
    print (radar_X.shape)

    opt_X = np.vstack((opt_X1,opt_X2))
    print (opt_X.shape)

    indices_X = np.vstack((indices_X1,indices_X2))
    print (indices_X.shape)

    y = np.vstack((y1,y2))
    print (y.shape)
    y = y[:,-1]

    print ('RADAR')
    bands = ["VH","VV"]
    simple_lr(radar_X,y,37,bands)
    print ('OPTICAL')
    bands = ["B2","B3","B4","B8","B5","B6","B7","B8A","B11","B12"]
    simple_lr(opt_X,y,37,bands)
    print ('INDICES')
    bands = ["NDVI","NDWI","EVI","MSAVI2","GDVI","CIGreen","CIRedEdge"]
    simple_lr(indices_X,y,37,bands)