R package

All the dependencies and needed functions to downscalled your data via SVM are contained inside an R package called “LETGmod”. It is up on github but not public, this why i have send you a line of code including a token to access it. To install it you have to install in the first place the package “devtools”, and then download it via install_github().

install.packages("devtools")
devtools::install_github("AlbanUR2/climathon@testing", auth_token = "SECRET_TOKEN")

Then you can just load the package as a normal one.

library(LETGmod)

Documentation

The package includes a partial documentation, some functions are associated with their documentation but the majority are not, and the vignette of the package is lacking. If you need to access the documentation of a function do so as normal in R.

?winkler_raster_HN

Files

The methodology is focused to the downscalling the minimum and maximun temperatures with the help of the topographic variables. The downscalling functions are designed to work with data.frame object to create the models and raster stack to assess the predictions. So in order to begging we will load our datas.

Extracting the topographic attributes

To extract the topographic attributes, you need a SRTM as a raster object and the extract_attributes() function.

srtm <- raster::raster("D:/test_application/MNT/mnt_25_l93.tif")
dem <- extract_attributes(dem = srtm)

Then you can extract the attributes value at each point of your data for that So you will need a spatial points for each oh your sampled individuals. We convert the extraction as data.frame for the following step.

variable <- raster::extract(dem, pts)
variable <- as.data.frame(variable)

Final files

In this tutorial my files are already created so i will load them into R:

tmax <- read.csv("D:/test_application/tmax_2012.csv", row.names = 1)
tmin <- read.csv("D:/test_application/tmin_2012.csv", row.names = 1)
variable <- read.csv("D:/test_application/variable_2012.csv", row.names = 1)
dem <- list.files("D:/test_application/MNT/", pattern = ".tif", full.names = TRUE)
dem <- raster::stack(dem)

Caution: The tmin and tmax colnames data.frame should be names like “XCapteur” for the names of your sensor is you have it, and “X%Year.%month.%%Day” for the following columns. Each columns correspond to the temperature for one day, rows are your sampled data.

The rows of data.frame “variable” must correspond to the rows of data.frame “tmin” and “tmax”. If you followed the method of extraction explained earlier you are ok.

Table output of tmin data.frame
XCapteur X2012.01.20 X2012.01.21 X2012.01.22 X2012.01.23 X2012.01.24
X1 9.196 9.571 10.199 8.199 3.437
X10 NA NA NA NA NA
X11 9.086 9.439 10.562 8.494 5.264
X12 NA NA NA NA NA
X14 8.801 9.239 9.923 7.998 5.019

Caution: Be sure that the variables names in the data.frame variable are matching the attributes names in raster.stack

Table output of variable data.frame
mnt_25_l93 pente_degree u v x y
34 0.9058503 -0.3162278 -0.9486833 6429775 448225
36 1.2809591 -0.8944272 0.4472136 6429150 448400
35 0.4051356 -0.7071068 -0.7071068 6429150 448550
38 0.0000000 0.0000000 1.0000000 6429400 448575
37 0.9058503 -0.9486833 -0.3162278 6429150 448700 
names(dem)
## [1] "mnt_25_l93"   "pente_degree" "u"            "v"            "x"           
## [6] "y"

Downscalling

SVM models

To create the downscalling model, just use the SVM_parral() function, you need to supply the temperatures of a day to downscall for Yvar argument, and the associated topographic attributes for each rows. Other arguments are aimed to parametization of the models for the e1071 package.

  model_day <-  SVM_parral(Yvar = tmin$X2012.01.01, variable = variable, k = 5, cost = 1, epsilon = 0.1, cross = 5, kernel = "radial")
## Time difference of 1.623566 secs

Predict

Using this model you can now predict the temperature for a whole area using the raster.stack of the topographic variables with the product_raster() function. The model’s object is contained in “model_day$svm_model”.

X2012.01.01_downscalled <- product_raster(model = model_day$svm_model, predictor = dem)
## Time difference of 5.546279 secs
raster::plot(X2012.01.01_downscalled)

Vegetative year workflow

The process to downscall a whole vegetative year is quite straight forward. You need to loop to downscalled tmin and tmax and for each one you can use lapply to downscalled each day. Here an example starting with models:

temperatures <- list()
temperatures[[1]] <- tmin[, -1]
temperatures[[2]] <- tmax[, -1]

models <- list()
for (i in 1:2) {
  models[[i]] <- lapply(temperatures[[i]], fun = function(x){
    output = SVM_parral(Yvar = x, variable = variable, k = 5, cost = 1, epsilon = 0.1, cross = 5, kernel = "radial")
    })
names(models[[i]]) <- colnames(temperatures[[i]])
  }

Then the predictions:

temperatures_downscalled <- list()

for (i in 1:2) {
  
  mode_svm <- list()
  
  for (j in 1:length(models[[i]])){
    mode_svm[[j]] <- models[[i]][[j]]$svm_model
  }
  names(mode_svm) <- names(models)
      
      
  temperatures_downscalled[[i]] = lapply(mode_svm, fun = function(y){ 
    prod <- product_raster(model = y, predictor = dem)
    return(prod)
      })
}

Index

Finally you can compute the bioclimatic indexes (i.e. Winkler, Huglin, GFV) using the output raster.stack with corresponding functions (i.e. winkler_raster_HN(), huglin_raster_HN(), Calcul_somme_GFV()).

The first two arguments of these functions are tmin and tmax stacks in that order, then the hemisphere (used to select the dates of the vegetative season) The huglin function have a k argument, choose it depending of lattitude of the data area. The GFV function have a Value_GFV argument, the critical degree-day sum paramter depending of the wanted cepage (see PARKER A.K. et al. 2011).

winkler <- winkler_raster_HN(temperatures_downscalled[[1]], temperatures_downscalled[[2]], hemisphere = "north")
huglin <- huglin_raster_HN(temperatures_downscalled[[1]], temperatures_downscalled[[2]], hemisphere = "north", k = 1.06)
gfv <- Calcul_somme_GFV(temperatures_downscalled[[1]], temperatures_downscalled[[2]], hemisphere = "north", Value_GFV = 2507)
mapview::mapview(winkler)
mapview::mapview(huglin)
# raster::plot(gfv)

You can also save the day by day sum of the winkler index to monitor the degree-day using the kee_means argument.

winkler <- winkler_raster_HN(temperatures_downscalled[[1]], temperatures_downscalled[[2]], hemisphere = "north", keep_means = TRUE)

winkler_raster_HN() and huglin_raster_HN() also include a method to compute the indexes for data.frame.

winkler_df <- winkler_raster_HN(tmin, tmax, hemisphere = "north")
huglin_df <- huglin_raster_HN(tmin, tmax, hemisphere = "north", k = 1.06)
winkler_df
##  [1] 1821.062 1743.185 1855.088       NA 1892.964 1744.310       NA 1712.335
##  [9] 1719.463 1740.496 1857.175       NA 1755.391 1659.987 1725.697       NA
## [17]       NA 1787.076 1717.425 1722.914       NA 1684.456 1732.457 1726.742
## [25]       NA       NA       NA 1811.003 1795.098 1777.785 1819.268       NA
## [33]       NA 1866.179       NA 1742.222 1872.108       NA 1733.958 1790.631
## [41] 1777.583 1745.812 1830.457 1875.314       NA 1845.210       NA 1755.090
## [49] 1748.338       NA 1718.686 1701.032       NA 1827.067       NA 1675.499
## [57] 1731.832 1690.054 1862.197       NA 1822.067 1702.423 1793.090 1732.110
## [65] 1626.829       NA       NA 1790.539       NA 1675.082 1798.898 1756.807
## [73] 1797.505 1792.482       NA 1680.284       NA 1871.438 1782.120 1739.780
## [81] 1749.618       NA       NA       NA       NA       NA       NA
huglin_df
##  [1] 2382.330 2297.987 2352.263 2322.720 2439.963 2318.312 2351.032 2343.831
##  [9] 2322.783 2332.901 2402.314 2243.576 2330.374 2286.713 2378.007 2431.392
## [17] 2404.469 2368.111 2320.880 2289.401       NA 2300.563 2377.762 2343.421
## [25] 2472.199 2363.105 2366.604 2374.577 2379.345 2388.099 2374.433 2337.022
## [33]       NA 2391.724 2300.216 2383.097 2403.793       NA 2248.824 2411.616
## [41] 2371.783 2316.321 2419.098 2416.616 2344.872 2408.984       NA 2282.930
## [49] 2322.172 2245.961 2269.454 2256.944 2378.188 2396.151       NA 2294.777
## [57] 2305.942 2294.987 2504.653 2337.048 2367.102 2300.055 2294.443 2311.488
## [65] 2279.560 2278.420 2295.746 2354.255 2370.499 2260.186 2405.854 2317.671
## [73] 2393.894 2357.419 2325.527 2295.209 2310.914 2434.317 2364.252 2298.073
## [81] 2323.875       NA       NA       NA       NA       NA       NA

Bonus Functions

Summarize tmin and tmax data.frame in data.frame ready to plot:

gg_temp <- summary_degree(tmin, tmax, "04.01", "09.30")
gg_temp$day <- as.Date(gg_temp$day, format = "X%Y.%m.%d")
ggplot(data = gg_temp, aes(y = mean, x = day, color = type)) +
  geom_line()

Transform data.frame who have day temperature separate in column, and a column for sensor name in standard data.format:

gg_min <- reformat_temp(tmin, "XCapteur")
ggplot(data = gg_min, aes(y = Temperature, x = Date, color = Capteur)) +
  geom_line()
## Warning: Removed 3149 row(s) containing missing values (geom_path).

Or conbine:

winkler_details <- winkler_raster_HN(tmin, tmax, hemisphere = "north", keep_means = TRUE)

degree_day <- as.data.frame(winkler_details$means[, 1])
for (i in 2:ncol(winkler_details$means)) {
  temp <- degree_day[, i-1] + winkler_details$means[, i]
  degree_day <- cbind(degree_day, temp)
}
colnames(degree_day) <- colnames(winkler_details$means)
degree_day$XCapteur <- tmin$XCapteur

gg_degre_day <- reformat_temp(degree_day, "XCapteur")
gg_degree_day_plty <- ggplot(data = gg_degre_day, aes(y = Temperature, x = Date, color = Capteur)) +
  geom_line()

ggplotly(gg_degree_day_plty)