Exploring regression-krigging models to map soil properties - Version 0.1
Ivan Lizarazo and Sebastian Gutierrez
6 de julio de 2018
Introduction
This is a preliminary report of a practical exercise exploring regression-krigging models for estimating soil properties from climatic variables.
Load libraries
## This script explores RK models from climate covariates
knitr::opts_chunk$set(fig.width = 8, collapse = TRUE)
setwd("C:\\Users\\UN-2344980\\Documents\\UN\\tallerDSM")
library(rgdal)## Loading required package: sp## rgdal: version: 1.3-3, (SVN revision 759)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/UN-2344980/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/UN-2344980/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1library(raster)
library(ranger)
library(GSIF)## GSIF version 0.5-4 (2017-04-25)## URL: http://gsif.r-forge.r-project.org/library(plotKML)## plotKML version 0.5-8 (2017-05-12)## URL: http://plotkml.r-forge.r-project.org/Study Area
Study area is the Sabana de Bogota, a montane savanna located in the southwestern part of the Altiplano Cundiboyacense in the center of Colombia.
library(ggplot2)
zona <- shapefile("C://Users//UN-2344980//Documents//UN//tallerDSM//gzona.shp")
zona
## class : SpatialPolygonsDataFrame
## features : 1
## extent : -74.56736, -72.66687, 4.015583, 6.117538 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 4
## names : Id, Area, BUFF_DIST, ORIG_FID
## value : 0, 1820674, 2000, 0
plot(zona, border="red", lwd=3, main = "Study Area")
zona_fortify <- fortify(zona)
## Regions defined for each Polygons
gg <- ggplot()
gg <- gg + geom_polygon(data=zona_fortify, aes(x=long, y=lat, group=group, fill=NA), color = "black", fill=NA, size=0.5)
gg <- gg <- gg + coord_map()
gg <- gg + labs(x=NULL, y=NULL,
title="Sabana de Bogota",
subtitle=NULL,
caption="Source: corposavia.gov.co")
gg
## Training samples
Training samples for this exercise were pre-processed by Viviana Varon and Gustavo Araujo from AGROSAVIA. Soil properties between 0-30 cm depth were calculated for each of the 267 sampled points. Such properties include Clase Textural, pH, CO, and CEC. A so-called Indice de Calidad de Suelos (SQI) was also computed using a weighted average of the above mentioned four variables.
load("C://Users//UN-2344980//Documents//UN//tallerDSM//raster_stack_001//matriz.R")
#summary(reg_matriz)
fort_matriz <- fortify(as.data.frame(reg_matriz))
zona_proj <- spTransform(zona,
crs(reg_matriz))
fzona <- fortify(zona_proj)
## Regions defined for each Polygons
coordinates(fzona) ~ long + lat
## coordinates(fzona) ~ long + lat
gg <- ggplot()
#gg <- gg + geom_polygon(data=fzona, aes(x=long, y=lat, group=group, fill=NA),
# color = "black", fill=NA, size=0.5)
gg <- gg + geom_point(data=fort_matriz, aes(x=XF, y=YF, color="red"))
#gg <- gg + coord_map()
gg <- gg + labs(x=NULL, y=NULL,
title="Sampled points",
subtitle=NULL,
caption="Source: corposavia.gov.co")
print(gg)
## Climatic covariates
WordClim-2 monthly precipitation, mean temperature and solar radiation will be tested as covariates for estimation of SQI, Texture, pH, CO and CEC.
precip <- list.files(path=".\\climate",pattern="*.tif",full.names=TRUE,recursive=TRUE)
precip
## [1] ".\\climate/Col_prec_01.tif" ".\\climate/Col_prec_02.tif"
## [3] ".\\climate/Col_prec_03.tif" ".\\climate/Col_prec_04.tif"
## [5] ".\\climate/Col_prec_05.tif" ".\\climate/Col_prec_06.tif"
## [7] ".\\climate/Col_prec_07.tif" ".\\climate/Col_prec_08.tif"
## [9] ".\\climate/Col_prec_09.tif" ".\\climate/Col_prec_10.tif"
## [11] ".\\climate/Col_prec_11.tif" ".\\climate/Col_prec_12.tif"
## [13] ".\\climate/col_srad_01.tif" ".\\climate/col_srad_02.tif"
## [15] ".\\climate/col_srad_03.tif" ".\\climate/col_srad_04.tif"
## [17] ".\\climate/col_srad_05.tif" ".\\climate/col_srad_06.tif"
## [19] ".\\climate/col_srad_07.tif" ".\\climate/col_srad_08.tif"
## [21] ".\\climate/col_srad_09.tif" ".\\climate/col_srad_10.tif"
## [23] ".\\climate/col_srad_11.tif" ".\\climate/col_srad_12.tif"
## [25] ".\\climate/tem_max_col_10.tif" ".\\climate/tem_max_col_2.tif"
## [27] ".\\climate/tem_max_col_4.tif" ".\\climate/tem_max_col_5.tif"
## [29] ".\\climate/tem_max_col_7.tif" ".\\climate/temp_max_col_1.tif"
## [31] ".\\climate/temp_max_col_11.tif" ".\\climate/temp_max_col_12.tif"
## [33] ".\\climate/temp_max_col_3.tif" ".\\climate/temp_max_col_6.tif"
## [35] ".\\climate/temp_max_col_8.tif" ".\\climate/temp_max_col_9.tif"
#env_data<- stack(results)
sprecip <- stack(precip) #stack all data
sprecip
## class : RasterStack
## dimensions : 251, 227, 56977, 36 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -74.56618, -72.67451, 4.018677, 6.110343 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : Col_prec_01, Col_prec_02, Col_prec_03, Col_prec_04, Col_prec_05, Col_prec_06, Col_prec_07, Col_prec_08, Col_prec_09, Col_prec_10, Col_prec_11, Col_prec_12, col_srad_01, col_srad_02, col_srad_03, ...
## min values : 1.0, 10.0, 39.0, 96.0, 80.0, 67.0, 25.0, 38.0, 54.0, 81.0, 44.0, 8.0, 15308.0, 14844.0, 15574.0, ...
## max values : 186.0, 231.0, 274.0, 469.0, 631.0, 745.0, 611.0, 554.0, 478.0, 489.0, 422.0, 335.0, 19080.0, 20880.0, 20284.0, ...
sprecip1 <- sprecip[[1]]
names(sprecip)
## [1] "Col_prec_01" "Col_prec_02" "Col_prec_03"
## [4] "Col_prec_04" "Col_prec_05" "Col_prec_06"
## [7] "Col_prec_07" "Col_prec_08" "Col_prec_09"
## [10] "Col_prec_10" "Col_prec_11" "Col_prec_12"
## [13] "col_srad_01" "col_srad_02" "col_srad_03"
## [16] "col_srad_04" "col_srad_05" "col_srad_06"
## [19] "col_srad_07" "col_srad_08" "col_srad_09"
## [22] "col_srad_10" "col_srad_11" "col_srad_12"
## [25] "tem_max_col_10" "tem_max_col_2" "tem_max_col_4"
## [28] "tem_max_col_5" "tem_max_col_7" "temp_max_col_1"
## [31] "temp_max_col_11" "temp_max_col_12" "temp_max_col_3"
## [34] "temp_max_col_6" "temp_max_col_8" "temp_max_col_9"Please note that the SpatialPointDataFrame which stores sampling points are in plane coordinates. On the contrary, the climatic covariates are in geographic coordinates.
Reprojecting and resampling covariates
Some pre-processing work is needed.
dem250 <- raster("C:\\Users\\UN-2344980\\Documents\\UN\\tallerDSM\\cov_topo\\dem_250.tif")
dem250
## class : RasterLayer
## dimensions : 932, 846, 788472 (nrow, ncol, ncell)
## resolution : 250, 250 (x, y)
## extent : 945443.1, 1156943, 935649.1, 1168649 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : C:\Users\UN-2344980\Documents\UN\tallerDSM\COV_TOPO\dem_250.tif
## names : dem_250
## values : 421, 4228 (min, max)
newcrs <- crs(dem250)
sprecip_proj <- projectRaster(sprecip, crs=newcrs) ## res=250, filename="npp_250.tif")
sprecip_proj
## class : RasterBrick
## dimensions : 261, 237, 61857, 36 (nrow, ncol, ncell, nlayers)
## resolution : 924, 922 (x, y)
## extent : 941113.7, 1160102, 931610.4, 1172252 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : in memory
## names : Col_prec_01, Col_prec_02, Col_prec_03, Col_prec_04, Col_prec_05, Col_prec_06, Col_prec_07, Col_prec_08, Col_prec_09, Col_prec_10, Col_prec_11, Col_prec_12, col_srad_01, col_srad_02, col_srad_03, ...
## min values : 1.203050, 10.000560, 39.033080, 96.435725, 81.000000, 67.906184, 27.251568, 36.277891, 54.137616, 81.214110, 46.030537, 8.299183, 15337.227914, 14981.272956, 15574.382807, ...
## max values : 204.82098, 222.22037, 268.12363, 468.05888, 615.46608, 727.53374, 601.45405, 550.49005, 466.82308, 489.07941, 416.25002, 331.21964, 19058.15522, 20867.82822, 20286.13104, ...
sprecip250 <- resample(sprecip_proj, dem250, method="ngb", filename="ncov_250.tif")
sprecip250
## class : RasterBrick
## dimensions : 932, 846, 788472, 36 (nrow, ncol, ncell, nlayers)
## resolution : 250, 250 (x, y)
## extent : 945443.1, 1156943, 935649.1, 1168649 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : in memory
## names : Col_prec_01, Col_prec_02, Col_prec_03, Col_prec_04, Col_prec_05, Col_prec_06, Col_prec_07, Col_prec_08, Col_prec_09, Col_prec_10, Col_prec_11, Col_prec_12, col_srad_01, col_srad_02, col_srad_03, ...
## min values : 1.203050, 10.000560, 39.033080, 96.435725, 81.000000, 67.906184, 27.251568, 36.277891, 54.137616, 81.214110, 46.030537, 8.299183, 15337.227914, 14981.272956, 15574.382807, ...
## max values : 204.82098, 222.22037, 268.12363, 468.05888, 615.46608, 727.53374, 601.45405, 550.49005, 466.82308, 489.07941, 416.25002, 331.21964, 19058.15522, 20867.82822, 20286.13104, ...
plot(sprecip250)
### Preparing regression matrix
Now, it’s time of preparing the regression matrix.
dim(reg_matriz)
## [1] 267 90
spoints <- reg_matriz[ -c(6:90) ]
#summary(spoints)
### extract precipitation values to sample points
dat <- extract(sprecip250, spoints, sp=TRUE)
#summary(dat)
mmatrix <- as.data.frame(dat)
dim(mmatrix)
## [1] 267 43
summary(mmatrix)
## Index_Cal Clase_textural pH CO
## Min. : 4.00 Lo :65 Min. :3.597 Min. : 0.060
## 1st Qu.:20.00 SaLo :52 1st Qu.:4.730 1st Qu.: 2.369
## Median :24.00 ClLo :50 Median :5.200 Median : 4.137
## Mean :24.43 Cl :41 Mean :5.246 Mean : 6.312
## 3rd Qu.:30.00 SaClLo :18 3rd Qu.:5.614 3rd Qu.: 7.928
## Max. :40.00 SiLo :13 Max. :7.726 Max. :40.472
## (Other):28
## CICA Col_prec_01 Col_prec_02 Col_prec_03
## Min. : 2.318 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.: 17.655 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median : 29.271 Median : 38.09 Median : 58.53 Median : 85.27
## Mean : 33.354 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.: 46.175 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :113.358 Max. :102.12 Max. :141.25 Max. :166.50
##
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
##
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
##
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
##
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
##
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
##
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
##
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
##
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
##
## temp_max_col_9 XF YF
## Min. : 8.892 Min. : 950168 Min. : 951586
## 1st Qu.:16.285 1st Qu.: 986051 1st Qu.:1014585
## Median :18.386 Median :1006549 Median :1034149
## Mean :17.666 Mean :1020551 Mean :1042184
## 3rd Qu.:19.285 3rd Qu.:1038248 3rd Qu.:1074343
## Max. :28.853 Max. :1146203 Max. :1164873
## Time for regression
Let’s create a regression model for estimation of the soil quality index (SQI). The first task is to refine the regression matrix:
newmatrix <- mmatrix[ -c(2:5) ]
dim(newmatrix)
## [1] 267 39
newmatrix <- newmatrix[ -c(38,39) ]
summary(newmatrix)
## Index_Cal Col_prec_01 Col_prec_02 Col_prec_03
## Min. : 4.00 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.:20.00 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median :24.00 Median : 38.09 Median : 58.53 Median : 85.27
## Mean :24.43 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.:30.00 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :40.00 Max. :102.12 Max. :141.25 Max. :166.50
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
## temp_max_col_9
## Min. : 8.892
## 1st Qu.:16.285
## Median :18.386
## Mean :17.666
## 3rd Qu.:19.285
## Max. :28.853
newmatrix <- newmatrix[complete.cases(newmatrix),]The second task is to create a Random Forest regression model to estimate the SQI:
library(ranger)
m.SQI <- ranger(Index_Cal ~ ., data= newmatrix, quantreg=TRUE, num.trees=150, seed=1234,
importance= 'impurity')
m.SQI
## Ranger result
##
## Call:
## ranger(Index_Cal ~ ., data = newmatrix, quantreg = TRUE, num.trees = 150, seed = 1234, importance = "impurity")
##
## Type: Regression
## Number of trees: 150
## Sample size: 267
## Number of independent variables: 36
## Mtry: 6
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 44.79973
## R squared (OOB): 0.0364456Based on the R squared value, it seems that this model is useless. Let’s try now a regression-krigging technique using the GSIF library:
mmatrix <- as.data.frame(dat)
dim(mmatrix)
## [1] 267 43
#summary(mmatrix)
newmatrix <- mmatrix[ -c(2:5)]
dim(newmatrix)
## [1] 267 39
#summary(newmatrix)
names(newmatrix)[38] <- "x"
names(newmatrix)[39] <- "y"
names(newmatrix)
## [1] "Index_Cal" "Col_prec_01" "Col_prec_02"
## [4] "Col_prec_03" "Col_prec_04" "Col_prec_05"
## [7] "Col_prec_06" "Col_prec_07" "Col_prec_08"
## [10] "Col_prec_09" "Col_prec_10" "Col_prec_11"
## [13] "Col_prec_12" "col_srad_01" "col_srad_02"
## [16] "col_srad_03" "col_srad_04" "col_srad_05"
## [19] "col_srad_06" "col_srad_07" "col_srad_08"
## [22] "col_srad_09" "col_srad_10" "col_srad_11"
## [25] "col_srad_12" "tem_max_col_10" "tem_max_col_2"
## [28] "tem_max_col_4" "tem_max_col_5" "tem_max_col_7"
## [31] "temp_max_col_1" "temp_max_col_11" "temp_max_col_12"
## [34] "temp_max_col_3" "temp_max_col_6" "temp_max_col_8"
## [37] "temp_max_col_9" "x" "y"
#newmatrix <- newmatrix[ -c(36,37) ]
#summary(newmatrix)
newmatrix <- newmatrix[complete.cases(newmatrix),]
#summary(newmatrix)
#sprecip250
pd <- as(sprecip250, "SpatialPixelsDataFrame")
#pd
plot(pd[2])
library(GSIF)
m2.SQI <- fit.regModel(Index_Cal ~ Col_prec_01 + Col_prec_02+ Col_prec_03+ Col_prec_04+ Col_prec_05+ Col_prec_06 + Col_prec_07
+ Col_prec_08 + Col_prec_09 + Col_prec_10+ Col_prec_11+ Col_prec_12 +col_srad_01 + col_srad_02 + col_srad_03 + col_srad_04 +
col_srad_05 + col_srad_06 +col_srad_07 +
col_srad_08 + col_srad_09 + col_srad_10 + col_srad_11 + col_srad_12
+ temp_max_col_1 + tem_max_col_2 + temp_max_col_3 +
tem_max_col_4 + tem_max_col_5 + temp_max_col_6
+ tem_max_col_7 + temp_max_col_8 + temp_max_col_9
+ tem_max_col_10 + temp_max_col_11 + temp_max_col_12 , rmatrix= newmatrix, pd,
method = "ranger")
## Fitting a randomForest model...
## Warning: Shapiro-Wilk normality test and Anderson-Darling normality test
## report probability of < .05 indicating lack of normal distribution for
## residuals
## Fitting a 2D variogram...
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Saving an object of class 'gstatModel'...
plot(m2.SQI)
#m2.SQIHow good is the regression-krigging model? Let’see.
(rk2 <- predict(m2.SQI, pd, nfold = 5))
## Generating predictions using the trend model (RK method)...
## Creating an object of class "SpatialPredictions"
## Variable : Index_Cal
## Minium value : 4
## Maximum value : 40
## Size : 267
## Total area : 48557125000
## Total area (units) : square-m
## Resolution (x) : 250
## Resolution (y) : 250
## Resolution (units) : m
## Vgm model : Exp
## Nugget (residual) : 14.4
## Sill (residual) : 0
## Range (residual) : 32200
## RMSE (validation) : NA
## Var explained : NA%
## Effective bytes : 62
## Compression method : gzipLet´s predict a SQI surface and evaluate observed vs predicted SQI values:
#(rk2 <- predict(m2.SQI, pd, nfold = 5))
slotNames(rk2)
## [1] "variable" "observed" "regModel.summary"
## [4] "vgmModel" "predicted" "validation"
e <- extract(stack(rk2@predicted), rk2@observed)
dim(e)
## [1] 267 2
names(data.frame(e))
## [1] "var1.var" "Index_Cal"
df <- cbind(rk2@observed@data, data.frame(e))
names(df)
## [1] "Index_Cal" "Col_prec_01" "Col_prec_02"
## [4] "Col_prec_03" "Col_prec_04" "Col_prec_05"
## [7] "Col_prec_06" "Col_prec_07" "Col_prec_08"
## [10] "Col_prec_09" "Col_prec_10" "Col_prec_11"
## [13] "Col_prec_12" "col_srad_01" "col_srad_02"
## [16] "col_srad_03" "col_srad_04" "col_srad_05"
## [19] "col_srad_06" "col_srad_07" "col_srad_08"
## [22] "col_srad_09" "col_srad_10" "col_srad_11"
## [25] "col_srad_12" "tem_max_col_10" "tem_max_col_2"
## [28] "tem_max_col_4" "tem_max_col_5" "tem_max_col_7"
## [31] "temp_max_col_1" "temp_max_col_11" "temp_max_col_12"
## [34] "temp_max_col_3" "temp_max_col_6" "temp_max_col_8"
## [37] "temp_max_col_9" "Index_Cal.residual" "Index_Cal.modelFit"
## [40] "var1.var" "Index_Cal"
cor(df[,1], df[,41])
## [1] 0.9131215
plotKML(rk2)
## KML file opened for writing...
## Reprojecting to +proj=longlat +datum=WGS84 ...
## Search for SAGA command line program and modules...
## Done
## Search for SAGA command line program and modules...
## Done
## Search for SAGA command line program and modules...
## Done
## Search for SAGA command line program and modules...
## Done
## Writing to KML...
## Reprojecting to +proj=longlat +datum=WGS84 ...
## Writing to KML...
## Closing rk2.kml
show(rk2)
## Variable : Index_Cal
## Minium value : 4
## Maximum value : 40
## Size : 267
## Total area : 48557125000
## Total area (units) : square-m
## Resolution (x) : 250
## Resolution (y) : 250
## Resolution (units) : m
## Vgm model : Exp
## Nugget (residual) : 14.4
## Sill (residual) : 0
## Range (residual) : 32200
## RMSE (validation) : NA
## Var explained : NA%
## Effective bytes : 62
## Compression method : gzip
summary(rk2)
## $variable
## [1] "Index_Cal"
##
## $minium
## [1] 4
##
## $maximum
## [1] 40
##
## $npoints
## [1] 267
##
## $area
## [1] "48557125000"
##
## $area.units
## [1] "square-m"
##
## $covariates
## [1] NA
##
## $family
## [1] NA
##
## $RMSE
## [1] NA
##
## $tvar
## [1] NA
##
## $breaks
## [1] NA
##
## $bonds
## [1] NA
##
## $Bytes
## [1] NA
##
## $compress
## [1] NA
##
## $npixels
## [1] 776914random Forest-based pH model
Let’s create several regression models for estimation of soil pH. Again, the first task is to refine the regression matrix. Next task is the creation of the model:
dim(mmatrix)
## [1] 267 43
newmatrix <- mmatrix[ -c(1,2,4,5) ]
dim(newmatrix)
## [1] 267 39
##summary(newmatrix)
newmatrix <- newmatrix[ -c(38,39) ]
summary(newmatrix)
## pH Col_prec_01 Col_prec_02 Col_prec_03
## Min. :3.597 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.:4.730 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median :5.200 Median : 38.09 Median : 58.53 Median : 85.27
## Mean :5.246 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.:5.614 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :7.726 Max. :102.12 Max. :141.25 Max. :166.50
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
## temp_max_col_9
## Min. : 8.892
## 1st Qu.:16.285
## Median :18.386
## Mean :17.666
## 3rd Qu.:19.285
## Max. :28.853
newmatrix <- newmatrix[complete.cases(newmatrix),]
m.pH <- ranger(pH ~ ., data= newmatrix, quantreg=TRUE, num.trees=150, seed=1234,
importance= 'impurity')
m.pH
## Ranger result
##
## Call:
## ranger(pH ~ ., data = newmatrix, quantreg = TRUE, num.trees = 150, seed = 1234, importance = "impurity")
##
## Type: Regression
## Number of trees: 150
## Sample size: 267
## Number of independent variables: 36
## Mtry: 6
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 0.389546
## R squared (OOB): 0.1924343random Forest-based CO model
Now, the turn is for CO estimation.
dim(mmatrix)
## [1] 267 43
newmatrix <- mmatrix[ -c(1,2,3,5) ]
dim(newmatrix)
## [1] 267 39
summary(newmatrix)
## CO Col_prec_01 Col_prec_02 Col_prec_03
## Min. : 0.060 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.: 2.369 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median : 4.137 Median : 38.09 Median : 58.53 Median : 85.27
## Mean : 6.312 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.: 7.928 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :40.472 Max. :102.12 Max. :141.25 Max. :166.50
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
## temp_max_col_9 XF YF
## Min. : 8.892 Min. : 950168 Min. : 951586
## 1st Qu.:16.285 1st Qu.: 986051 1st Qu.:1014585
## Median :18.386 Median :1006549 Median :1034149
## Mean :17.666 Mean :1020551 Mean :1042184
## 3rd Qu.:19.285 3rd Qu.:1038248 3rd Qu.:1074343
## Max. :28.853 Max. :1146203 Max. :1164873
newmatrix <- newmatrix[ -c(38,39) ]
summary(newmatrix)
## CO Col_prec_01 Col_prec_02 Col_prec_03
## Min. : 0.060 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.: 2.369 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median : 4.137 Median : 38.09 Median : 58.53 Median : 85.27
## Mean : 6.312 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.: 7.928 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :40.472 Max. :102.12 Max. :141.25 Max. :166.50
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
## temp_max_col_9
## Min. : 8.892
## 1st Qu.:16.285
## Median :18.386
## Mean :17.666
## 3rd Qu.:19.285
## Max. :28.853
newmatrix <- newmatrix[complete.cases(newmatrix),]
m.CO <- ranger(CO ~ ., data= newmatrix, quantreg=TRUE, num.trees=150, seed=1234,
importance= 'impurity')
m.CO
## Ranger result
##
## Call:
## ranger(CO ~ ., data = newmatrix, quantreg = TRUE, num.trees = 150, seed = 1234, importance = "impurity")
##
## Type: Regression
## Number of trees: 150
## Sample size: 267
## Number of independent variables: 36
## Mtry: 6
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 24.42829
## R squared (OOB): 0.3283947a regression-kriging CO model
dim(mmatrix)
## [1] 267 43
names(mmatrix)
## [1] "Index_Cal" "Clase_textural" "pH"
## [4] "CO" "CICA" "Col_prec_01"
## [7] "Col_prec_02" "Col_prec_03" "Col_prec_04"
## [10] "Col_prec_05" "Col_prec_06" "Col_prec_07"
## [13] "Col_prec_08" "Col_prec_09" "Col_prec_10"
## [16] "Col_prec_11" "Col_prec_12" "col_srad_01"
## [19] "col_srad_02" "col_srad_03" "col_srad_04"
## [22] "col_srad_05" "col_srad_06" "col_srad_07"
## [25] "col_srad_08" "col_srad_09" "col_srad_10"
## [28] "col_srad_11" "col_srad_12" "tem_max_col_10"
## [31] "tem_max_col_2" "tem_max_col_4" "tem_max_col_5"
## [34] "tem_max_col_7" "temp_max_col_1" "temp_max_col_11"
## [37] "temp_max_col_12" "temp_max_col_3" "temp_max_col_6"
## [40] "temp_max_col_8" "temp_max_col_9" "XF"
## [43] "YF"
newmatrix <- mmatrix[ -c(1,2,3,5) ]
dim(newmatrix)
## [1] 267 39
summary(newmatrix)
## CO Col_prec_01 Col_prec_02 Col_prec_03
## Min. : 0.060 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.: 2.369 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median : 4.137 Median : 38.09 Median : 58.53 Median : 85.27
## Mean : 6.312 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.: 7.928 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :40.472 Max. :102.12 Max. :141.25 Max. :166.50
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
## temp_max_col_9 XF YF
## Min. : 8.892 Min. : 950168 Min. : 951586
## 1st Qu.:16.285 1st Qu.: 986051 1st Qu.:1014585
## Median :18.386 Median :1006549 Median :1034149
## Mean :17.666 Mean :1020551 Mean :1042184
## 3rd Qu.:19.285 3rd Qu.:1038248 3rd Qu.:1074343
## Max. :28.853 Max. :1146203 Max. :1164873
names(newmatrix)[38] <- "x"
names(newmatrix)[39] <- "y"
names(newmatrix)
## [1] "CO" "Col_prec_01" "Col_prec_02"
## [4] "Col_prec_03" "Col_prec_04" "Col_prec_05"
## [7] "Col_prec_06" "Col_prec_07" "Col_prec_08"
## [10] "Col_prec_09" "Col_prec_10" "Col_prec_11"
## [13] "Col_prec_12" "col_srad_01" "col_srad_02"
## [16] "col_srad_03" "col_srad_04" "col_srad_05"
## [19] "col_srad_06" "col_srad_07" "col_srad_08"
## [22] "col_srad_09" "col_srad_10" "col_srad_11"
## [25] "col_srad_12" "tem_max_col_10" "tem_max_col_2"
## [28] "tem_max_col_4" "tem_max_col_5" "tem_max_col_7"
## [31] "temp_max_col_1" "temp_max_col_11" "temp_max_col_12"
## [34] "temp_max_col_3" "temp_max_col_6" "temp_max_col_8"
## [37] "temp_max_col_9" "x" "y"
#newmatrix <- newmatrix[ -c(36,37) ]
#summary(newmatrix)
newmatrix <- newmatrix[complete.cases(newmatrix),]
summary(newmatrix)
## CO Col_prec_01 Col_prec_02 Col_prec_03
## Min. : 0.060 Min. : 8.68 Min. : 15.37 Min. : 41.18
## 1st Qu.: 2.369 1st Qu.: 28.78 1st Qu.: 49.93 1st Qu.: 77.35
## Median : 4.137 Median : 38.09 Median : 58.53 Median : 85.27
## Mean : 6.312 Mean : 40.58 Mean : 62.24 Mean : 93.03
## 3rd Qu.: 7.928 3rd Qu.: 48.77 3rd Qu.: 72.04 3rd Qu.:108.46
## Max. :40.472 Max. :102.12 Max. :141.25 Max. :166.50
## Col_prec_04 Col_prec_05 Col_prec_06 Col_prec_07
## Min. :102.2 Min. : 90.86 Min. : 70.77 Min. : 35.99
## 1st Qu.:122.3 1st Qu.:118.50 1st Qu.:100.88 1st Qu.: 64.11
## Median :136.0 Median :136.86 Median :117.19 Median : 81.73
## Mean :143.3 Mean :142.53 Mean :137.91 Mean :107.26
## 3rd Qu.:157.1 3rd Qu.:162.53 3rd Qu.:153.57 3rd Qu.:107.81
## Max. :246.5 Max. :266.30 Max. :487.72 Max. :457.34
## Col_prec_08 Col_prec_09 Col_prec_10 Col_prec_11
## Min. : 45.00 Min. : 59.44 Min. : 87.07 Min. : 72.55
## 1st Qu.: 69.39 1st Qu.: 86.60 1st Qu.:141.22 1st Qu.:104.44
## Median : 80.47 Median :100.36 Median :154.85 Median :123.50
## Mean : 95.26 Mean :109.20 Mean :161.20 Mean :129.29
## 3rd Qu.:108.06 3rd Qu.:118.33 3rd Qu.:173.97 3rd Qu.:141.75
## Max. :289.89 Max. :238.61 Max. :313.68 Max. :293.40
## Col_prec_12 col_srad_01 col_srad_02 col_srad_03
## Min. : 22.97 Min. :16016 Min. :15449 Min. :16089
## 1st Qu.: 53.97 1st Qu.:17303 1st Qu.:17245 1st Qu.:17175
## Median : 71.12 Median :17532 Median :17912 Median :17568
## Mean : 73.04 Mean :17478 Mean :17934 Mean :17487
## 3rd Qu.: 92.53 3rd Qu.:17747 3rd Qu.:18578 3rd Qu.:17903
## Max. :155.82 Max. :18311 Max. :19767 Max. :19748
## col_srad_04 col_srad_05 col_srad_06 col_srad_07
## Min. :13894 Min. :14730 Min. :14171 Min. :14502
## 1st Qu.:15501 1st Qu.:15376 1st Qu.:15424 1st Qu.:15514
## Median :15795 Median :15475 Median :15586 Median :15983
## Mean :15594 Mean :15490 Mean :15569 Mean :15818
## 3rd Qu.:15890 3rd Qu.:15591 3rd Qu.:15737 3rd Qu.:16123
## Max. :17440 Max. :16333 Max. :17158 Max. :17112
## col_srad_08 col_srad_09 col_srad_10 col_srad_11
## Min. :15984 Min. :16480 Min. :14869 Min. :14594
## 1st Qu.:16580 1st Qu.:16870 1st Qu.:15902 1st Qu.:15464
## Median :16678 Median :17225 Median :16025 Median :15838
## Mean :16725 Mean :17199 Mean :15963 Mean :15723
## 3rd Qu.:16796 3rd Qu.:17473 3rd Qu.:16179 3rd Qu.:15991
## Max. :17854 Max. :17936 Max. :16364 Max. :16241
## col_srad_12 tem_max_col_10 tem_max_col_2 tem_max_col_4
## Min. :14688 Min. : 9.186 Min. : 8.385 Min. : 9.085
## 1st Qu.:15239 1st Qu.:16.277 1st Qu.:15.748 1st Qu.:16.126
## Median :15913 Median :18.300 Median :18.162 Median :18.453
## Mean :15954 Mean :17.662 Mean :17.277 Mean :17.579
## 3rd Qu.:16532 3rd Qu.:19.299 3rd Qu.:18.696 3rd Qu.:19.092
## Max. :17741 Max. :28.251 Max. :29.140 Max. :28.842
## tem_max_col_5 tem_max_col_7 temp_max_col_1 temp_max_col_11
## Min. : 8.686 Min. : 8.193 Min. : 8.398 Min. : 9.492
## 1st Qu.:15.786 1st Qu.:15.840 1st Qu.:15.481 1st Qu.:16.437
## Median :18.268 Median :18.193 Median :17.764 Median :18.377
## Mean :17.343 Mean :17.281 Mean :16.952 Mean :17.772
## 3rd Qu.:18.892 3rd Qu.:18.974 3rd Qu.:18.300 3rd Qu.:19.384
## Max. :28.740 Max. :28.642 Max. :28.951 Max. :28.051
## temp_max_col_12 temp_max_col_3 temp_max_col_6 temp_max_col_8
## Min. : 9.292 Min. : 8.986 Min. : 7.985 Min. : 8.786
## 1st Qu.:16.287 1st Qu.:16.276 1st Qu.:15.276 1st Qu.:16.300
## Median :18.278 Median :18.540 Median :17.889 Median :18.205
## Mean :17.628 Mean :17.685 Mean :16.816 Mean :17.520
## 3rd Qu.:19.195 3rd Qu.:19.160 3rd Qu.:18.405 3rd Qu.:19.080
## Max. :28.251 Max. :29.240 Max. :28.551 Max. :28.951
## temp_max_col_9 x y
## Min. : 8.892 Min. : 950168 Min. : 951586
## 1st Qu.:16.285 1st Qu.: 986051 1st Qu.:1014585
## Median :18.386 Median :1006549 Median :1034149
## Mean :17.666 Mean :1020551 Mean :1042184
## 3rd Qu.:19.285 3rd Qu.:1038248 3rd Qu.:1074343
## Max. :28.853 Max. :1146203 Max. :1164873
#sprecip250
pd <- as(sprecip250, "SpatialPixelsDataFrame")
#pdNow, the r-k model:
colnames(newmatrix)
## [1] "CO" "Col_prec_01" "Col_prec_02"
## [4] "Col_prec_03" "Col_prec_04" "Col_prec_05"
## [7] "Col_prec_06" "Col_prec_07" "Col_prec_08"
## [10] "Col_prec_09" "Col_prec_10" "Col_prec_11"
## [13] "Col_prec_12" "col_srad_01" "col_srad_02"
## [16] "col_srad_03" "col_srad_04" "col_srad_05"
## [19] "col_srad_06" "col_srad_07" "col_srad_08"
## [22] "col_srad_09" "col_srad_10" "col_srad_11"
## [25] "col_srad_12" "tem_max_col_10" "tem_max_col_2"
## [28] "tem_max_col_4" "tem_max_col_5" "tem_max_col_7"
## [31] "temp_max_col_1" "temp_max_col_11" "temp_max_col_12"
## [34] "temp_max_col_3" "temp_max_col_6" "temp_max_col_8"
## [37] "temp_max_col_9" "x" "y"
m3.CO <- fit.regModel(CO ~ Col_prec_01 + Col_prec_02+ Col_prec_03+ Col_prec_04+ Col_prec_05+ Col_prec_06 + Col_prec_07
+ Col_prec_08 + Col_prec_09 + Col_prec_10+ Col_prec_11+ Col_prec_12 +col_srad_01 + col_srad_02 + col_srad_03 + col_srad_04 +
col_srad_05 + col_srad_06 +col_srad_07 +
col_srad_08 + col_srad_09 + col_srad_10 + col_srad_11 + col_srad_12
+ temp_max_col_1 + tem_max_col_2 + temp_max_col_3 +
tem_max_col_4 + tem_max_col_5 + temp_max_col_6
+ tem_max_col_7 + temp_max_col_8 + temp_max_col_9
+ tem_max_col_10 + temp_max_col_11 + temp_max_col_12, rmatrix= newmatrix, pd,
method = "ranger")
## Fitting a randomForest model...
## Warning: Shapiro-Wilk normality test and Anderson-Darling normality test
## report probability of < .05 indicating lack of normal distribution for
## residuals
## Fitting a 2D variogram...
## Warning in gstat::fit.variogram(svgm, model = ivgm, ...): No convergence
## after 200 iterations: try different initial values?
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Saving an object of class 'gstatModel'...
plot(m3.CO)
(rk6 <- predict(m3.CO, pd, nfold = 5))
## Generating predictions using the trend model (RK method)...
## Creating an object of class "SpatialPredictions"
## Variable : CO
## Minium value : 0.06
## Maximum value : 40.47249
## Size : 267
## Total area : 48557125000
## Total area (units) : square-m
## Resolution (x) : 250
## Resolution (y) : 250
## Resolution (units) : m
## Vgm model : Exp
## Nugget (residual) : 4.45
## Sill (residual) : 0
## Range (residual) : 32200
## RMSE (validation) : NA
## Var explained : NA%
## Effective bytes : 62
## Compression method : gzip
slotNames(rk6)
## [1] "variable" "observed" "regModel.summary"
## [4] "vgmModel" "predicted" "validation"
e <- extract(stack(rk6@predicted), rk6@observed)
#e
names(data.frame(e))
## [1] "var1.var" "CO"
df <- cbind(rk6@observed@data, data.frame(e))
names(df)
## [1] "CO" "Col_prec_01" "Col_prec_02"
## [4] "Col_prec_03" "Col_prec_04" "Col_prec_05"
## [7] "Col_prec_06" "Col_prec_07" "Col_prec_08"
## [10] "Col_prec_09" "Col_prec_10" "Col_prec_11"
## [13] "Col_prec_12" "col_srad_01" "col_srad_02"
## [16] "col_srad_03" "col_srad_04" "col_srad_05"
## [19] "col_srad_06" "col_srad_07" "col_srad_08"
## [22] "col_srad_09" "col_srad_10" "col_srad_11"
## [25] "col_srad_12" "tem_max_col_10" "tem_max_col_2"
## [28] "tem_max_col_4" "tem_max_col_5" "tem_max_col_7"
## [31] "temp_max_col_1" "temp_max_col_11" "temp_max_col_12"
## [34] "temp_max_col_3" "temp_max_col_6" "temp_max_col_8"
## [37] "temp_max_col_9" "CO.residual" "CO.modelFit"
## [40] "var1.var" "CO"
cor(df[,1], df[,41])
## [1] 0.9440865
## [1] 0.9044326
plotKML(rk6)
## KML file opened for writing...
## Reprojecting to +proj=longlat +datum=WGS84 ...
## Search for SAGA command line program and modules...
## Done
## Search for SAGA command line program and modules...
## Done
## Search for SAGA command line program and modules...
## Done
## Search for SAGA command line program and modules...
## Done
## Writing to KML...
## Reprojecting to +proj=longlat +datum=WGS84 ...
## Writing to KML...
## Closing rk6.kml
show(rk6)
## Variable : CO
## Minium value : 0.06
## Maximum value : 40.47249
## Size : 267
## Total area : 48557125000
## Total area (units) : square-m
## Resolution (x) : 250
## Resolution (y) : 250
## Resolution (units) : m
## Vgm model : Exp
## Nugget (residual) : 4.45
## Sill (residual) : 0
## Range (residual) : 32200
## RMSE (validation) : NA
## Var explained : NA%
## Effective bytes : 62
## Compression method : gzip
summary(rk6)
## $variable
## [1] "CO"
##
## $minium
## [1] 0.06
##
## $maximum
## [1] 40.47249
##
## $npoints
## [1] 267
##
## $area
## [1] "48557125000"
##
## $area.units
## [1] "square-m"
##
## $covariates
## [1] NA
##
## $family
## [1] NA
##
## $RMSE
## [1] NA
##
## $tvar
## [1] NA
##
## $breaks
## [1] NA
##
## $bonds
## [1] NA
##
## $Bytes
## [1] NA
##
## $compress
## [1] NA
##
## $npixels
## [1] 776914
## writeRaster(rk6.kml, 'rk6.tif')
par(mfrow=c(1,2))
plot(pd)
plot(rk6)