Commit d8c36ec1 authored by luroth's avatar luroth
Browse files

begin stem elongation estimation

parent 6d12ce0b
......@@ -10,11 +10,14 @@ from scipy.ndimage.measurements import center_of_mass, label
from skimage.morphology import watershed
from scipy.sparse import csr_matrix
import math
from scipy.interpolate import UnivariateSpline
from scipy.interpolate import UnivariateSpline, interp1d
from sklearn.linear_model import LinearRegression
from matplotlib import pyplot as plt
from skimage.morphology import disk
from skimage.filters import rank
from sklearn import svm
from sklearn.preprocessing import StandardScaler
# Helper functions
......@@ -343,10 +346,16 @@ def process_NadirCC_LCCC_LA_AI(path_campaign, campaign_date, GSD):
# Trains SVM
df_training = pd.read_csv('./GroundAerialCoverage/svm_SE_training.csv')
features_col = [col for col in df_training if col.startswith('gt_bin')]
sc_x = StandardScaler()
X = df_training[features_col]
sc_x.fit(X)
X = sc_x.transform(X)
y = df_training['delta_to_BBCH']
svm_predictor = svm.SVR(C=150, kernel='rbf')
svm_predictor = svm.SVR(C=2, kernel='rbf', gamma = 0.01, epsilon = 0.1)
svm_predictor.fit(X, y)
print("Process campaign ", path_campaign)
......@@ -470,7 +479,7 @@ def process_NadirCC_LCCC_LA_AI(path_campaign, campaign_date, GSD):
X = [list(d.values()) for d in df_LA['gt_bins']]
value = -svm_predictor.predict(X)
value = -svm_predictor.predict(sc_x.transform(X))
df_LA['value'] = value
......@@ -480,7 +489,7 @@ def process_NadirCC_LCCC_LA_AI(path_campaign, campaign_date, GSD):
def process_begin_stem_elongation(paths_campaigns):
def process_begin_stem_elongation(paths_campaigns, graph=True):
......@@ -515,7 +524,7 @@ def process_begin_stem_elongation(paths_campaigns):
df_AI['value_abs'] = df_AI.value.abs()
# Find transition from negative to positive index (begin SE) with spline
# Find transition from negative to positive index (begin SE) with lm
df_AI_smooth = df_AI.copy()
df_AI_smooth['yday'] = pd.to_numeric(pd.to_datetime(df_AI_smooth.timestamp).dt.strftime('%j'), downcast="float")
df_AI_smooth['year'] = pd.to_numeric(pd.to_datetime(df_AI_smooth.timestamp).dt.strftime('%Y'))
......@@ -523,25 +532,46 @@ def process_begin_stem_elongation(paths_campaigns):
bbch30_dates = {}
plot_groups = df_AI_smooth.groupby("plot_label")
for plot_label, df_plot_group in plot_groups:
print("Smoothing plot", plot_label)
spline = UnivariateSpline(y=df_plot_group.value, x=df_plot_group.yday, w= -np.abs(df_plot_group.value), k=3)
roots = spline.roots()
first_root = spline.roots()[0] if len(roots) > 0 else np.max(df_plot_group.yday)
from matplotlib import pyplot as plt
plt.scatter(df_plot_group.yday, df_plot_group.value)
x_range = np.arange(np.min(df_plot_group.yday), np.max(df_plot_group.yday), 1)
plt.plot(x_range, spline(x_range), label='Cubic Spline')
plt.vlines(x=first_root, ymin=-1, ymax=1)
plt.title(plot_label)
#plt.show()
plt.savefig("/home/luroth/Downloads/figs/" + plot_label + ".png")
plt.close()
print("Estimating transition trough zero for plot", plot_label)
# For model: find closest point before zero and two closes points after zero
df_before_zero = df_plot_group[df_plot_group.value < -1].nlargest(1, 'campaign_date')
df_after_zero = df_plot_group[df_plot_group.value >= -1].nsmallest(2, 'campaign_date')
# If oldest after zero date is older than newest before zero date: add one before zero
i_before_zero = 1
while np.min(df_after_zero.campaign_date) < np.min(df_before_zero.campaign_date) and i_before_zero < 10:
i_before_zero += 1
df_before_zero = df_plot_group[(df_plot_group.value < -1) & (df_plot_group.campaign_date < np.min(df_before_zero.campaign_date))].nlargest(1, 'campaign_date')
df_plot_group_lm = pd.concat([df_before_zero, df_after_zero])
lm = None
if len(df_plot_group_lm) < 2:
# Set date to max date if less than two points for lm
lm_zero = np.max(df_plot_group.yday)
else:
# Build lm, calculate zero transition
lm = LinearRegression().fit([[v] for v in df_plot_group_lm.yday], df_plot_group_lm.value)
lm_zero_ = - lm.intercept_ / lm.coef_[0]
# If zero transition is outside max date: correct extrapolation
lm_zero = lm_zero_ if (lm_zero_ < np.max(df_plot_group.yday) and lm.coef_[0] > 0) else np.max(df_plot_group.yday)
# Plot
if graph:
plt.scatter(df_plot_group.yday, df_plot_group.value)
x_range = np.arange(np.min(df_plot_group.yday), np.max(df_plot_group.yday), 1)
plt.scatter(df_plot_group_lm.yday, df_plot_group_lm.value)
if lm:
plt.plot(x_range, x_range * lm.coef_[0] + lm.intercept_, label='lm')
plt.plot(lm_zero, -1, marker = 'x', c='black')
plt.hlines(y=-1, xmin=np.min(x_range), xmax=np.max(x_range))
plt.title(plot_label)
#plt.show()
plt.savefig("/home/luroth/Downloads/figs/" + plot_label + ".png")
plt.close()
# Convert back to date
year = np.unique(df_plot_group.year)[0]
bbch30_date = datetime(year, 1, 1) + timedelta(np.round(first_root) - 1)
bbch30_date = datetime(year, 1, 1) + timedelta(np.round(lm_zero) - 1)
bbch30_dates[plot_label] = bbch30_date
......
......@@ -2,7 +2,7 @@ from GroundAerialCoverage import CanopyAnalysis
from pathlib import Path
path_p = Path("P:")
#path_p = Path("/home/luroth/public")
path_p = Path("/home/luroth/public")
if __name__ == "__main__":
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment