Skip to content

Commit

Permalink
Adds logs to testthat file
Browse files Browse the repository at this point in the history
  • Loading branch information
vwmaus committed Oct 19, 2023
1 parent 74f341f commit c8489d2
Showing 1 changed file with 41 additions and 17 deletions.
58 changes: 41 additions & 17 deletions tests/testthat/test-twdtw_classify.R
Original file line number Diff line number Diff line change
@@ -1,28 +1,28 @@
cat("Loading mgcv...\n")
library(mgcv)

# Read training samples
samples <- st_read(system.file("mato_grosso_brazil/samples.gpkg", package = "dtwSat"), quiet = TRUE)
cat("Reading training samples...\n")
samples <- st_read(system.file("mato_grosso_brazil/samples.gpkg", package = "dtwSat"))

# Satellite image time sereis files
cat("Geting satellite image time series files...\n")
tif_files <- dir(system.file("mato_grosso_brazil", package = "dtwSat"), pattern = "\\.tif$", full.names = TRUE)

# The acquisition date is in the file name are not the true acquisition date
# of each pixel. MOD13Q1 is a 16-day composite product, so the acquisition date
# is the first day of the 16-day period
cat("Geting acquisition dates from file names...\n")
acquisition_date <- as.Date(regmatches(tif_files, regexpr("[0-9]{8}", tif_files)), format = "%Y%m%d")

# Read the data as a stars object setting the time/date for each observation
# using along. This will prodcue a 4D array (data-cube) which will then be converted
# to a 3D array by spliting the 'band' dimension
cat("Reading data as stars object...\n")
dc <- read_stars(tif_files,
proxy = FALSE,
along = list(time = acquisition_date),
RasterIO = list(bands = 1:6))

cat("Setting band names...\n")
dc <- st_set_dimensions(dc, 3, c("EVI", "NDVI", "RED", "BLUE", "NIR", "MIR"))

cat("Converting 4D into a 3D datacube...\n")
dc <- split(dc, c("band"))

# Create a knn1-twdtw model
cat("Creating twdtw_knn1 model using GAM...\n")
system.time(
m <- twdtw_knn1(x = dc,
y = samples,
Expand All @@ -32,22 +32,24 @@ system.time(
time_weight = c(steepness = 0.1, midpoint = 50))
)

cat("Printing model...\n")
print(m)

# Visualize model patterns
cat("Visualizing model patterns...\n")
plot(m)

# Classify satellite image time series
cat("Classifying satellite image time series...\n")
system.time(lu <- predict(dc, model = m))

cat("Printing land use classification...\n")
print(lu)

# # Visualise land use classification
cat("Visualizing land use classification...\n")
ggplot() +
geom_stars(data = lu) +
theme_minimal()

# ### OTHER TESTS
cat("Testing model with resampling frequency of 60 and GAM...\n")
m <- twdtw_knn1(x = dc,
y = samples,
cycle_length = 'year',
Expand All @@ -56,8 +58,10 @@ m <- twdtw_knn1(x = dc,
smooth_fun = function(x, y) gam(y ~ s(x), data = data.frame(x = x, y = y)),
resampling_freq = 60)

cat("Visualizing resampled patterns...\n")
plot(m)

cat("Classifying satellite image time series with resampled patterns...\n")
system.time(
lu <- predict(dc,
model = m,
Expand All @@ -66,37 +70,57 @@ system.time(
time_weight = c(steepness = 0.1, midpoint = 50))
)

cat("Printing land use classification of resampled patterns...\n")
print(lu)

# Visualise land use classification
cat("Visualizing land use classification of resampled patterns...\n")
ggplot() +
geom_stars(data = lu) +
theme_minimal()

# Test model without samples reduction
cat("Testing model without sample reduction...\n")
m <- twdtw_knn1(x = dc,
y = samples,
cycle_length = 'year',
time_scale = 'day',
time_weight = c(steepness = 0.1, midpoint = 50),
smooth_fun = NULL)

cat("Printing model without sample reduction...\n")
print(m)

cat("Visualizing model patterns without sample reduction...\n")
plot(m)

cat("Visualizing model patterns without sample reduction for EVI and NDVI...\n")
plot(m, bands = c('EVI', 'NDVI'))

# Test custom smooth function -- mean value for each time in dc
cat("Testing custom smooth function using the average for each observation date...\n")
m <- twdtw_knn1(x = dc,
y = samples,
cycle_length = 'year',
time_scale = 'day',
time_weight = c(steepness = 0.1, midpoint = 50),
smooth_fun = function(x, y) lm(y ~ factor(x), data = data.frame(x = x, y = y)))

cat("Printing model with custom smooth...\n")
print(m)

cat("Visualizing model patterns with custom smooth...\n")
plot(m)

cat("Visualizing model patterns with custom smooth for EVI and NDVI...\n")
plot(m, bands = c('EVI', 'NDVI'))

cat("Classifying satellite image time series with custom smooth...\n")
system.time(lu <- predict(dc, model = m))

cat("Printing land use classification with custom smooth...\n")
print(lu)

cat("Visualizing land use classification with custom smooth...\n")
ggplot() +
geom_stars(data = lu) +
theme_minimal()

cat("Script completed.\n")

0 comments on commit c8489d2

Please sign in to comment.