Detecting change point in a time series I'm dealing with time series from satellite imagery, where I have a sudden change (drop), that I can see from the plot, but I need a statistical test to detect it. I already checked for stationarity using Dickey-Fuller, and this test doesn't detect the change. 
Basically the time series has a first part with a certain average, and the second part with a lower average, and I need to detect where (when) the change occurs.
The time series also have seasonality.
Thanks for any pointers
Edit: this is a sample of the data, that may or may not have a change point: https://pastebin.com/bVSsxVtR
 A: You can do this with mcp. First, let's get your data in an accessible format. The variable "days" is the number of days since the first record. I remove the NAs:
library(dplyr)
D = read.table("bVSsxVtR.txt") %>%
  rename(row = V1, date = V2, y = V3) %>%
  mutate(days = as.numeric(as.POSIXct(date) - as.POSIXct("2014-10-31")) / (24*3600)) %>%
  select(-row)

head(D)

  V1       date         y days
1  1 2014-10-31 -18.13354    0
2  2 2014-11-04 -17.71648    4
3  4 2014-11-07 -15.03215    7
4  5 2014-11-24 -15.66677   24
5  6 2014-12-01 -18.37760   31
6  8 2014-12-13 -16.23460   43

Now let's fit to the data a model with a change point and a model without a change point. We add an AR(1) to both, to capture some of the time-series autocorrelation, but you can choose the order to your liking and add an ar(N) term to segment two too, if the change involves a change in the AR coefficient(s). Note that the dates are not equally spaced so the AR is not accurate, though it may be good enough.
model = list(y ~ 1 + ar(1), ~ 1)
fit = mcp(model, D, par_x = "days")

model_null = list(y ~ 1 + ar(1))
fit_null = mcp(model_null, D, par_x = "days")

Here's plot(fit) + plot(fit_null), i.e., the model that was forced to identify a change point and the one without. The posterior (blue distribution on the right-hand plot) shows that if there is a change point, its location is quite well defined:

Now for the critical part: we use leave-one-out cross-validation to compare the predictive performance of the two models. 
fit$loo = loo(fit)
fit_null$loo = loo(fit_null)
loo::loo_compare(fit$loo, fit_null$loo)

       elpd_diff se_diff
model1  0.0       0.0   
model2 -4.4       3.3  

The change point model is preferred, but only slightly (small elpd_diff / elpd_se ratio). You can read more about modeling time-series, model comparison, etc. on the mcp website.
