Loading layers
shp_aru = st_read('../../arauca/81_ARAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp')
## Reading layer `MGN_MPIO_POLITICO' from data source
## `C:\Users\ADMIN\Documentos\CARLOS\UNAL\Asignaturas\Geomatica\arauca\81_ARAUCA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 10 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -72.36662 ymin: 6.036228 xmax: -69.42756 ymax: 7.104381
## Geodetic CRS: WGS 84
shp_aru_magna = st_transform(shp_aru, crs = 3116)
Descargar las 4 partes del raster de SoilGrids
# bbox = as.vector(st_bbox(shp_aru))
# mx = mean(bbox[c(1,3)])
# my = mean(bbox[c(2,4)])
#
# f_part2_url <- function(cord){
# paste0('long(',cord[1],',',cord[2],')&SUBSET=lat(',cord[3],',',cord[4],')')
# }
#
# ras_coor = rbind(
# c(bbox[1], mx, bbox[2], my),
# c(mx, bbox[3], bbox[2], my),
# c(bbox[1], mx, my, bbox[4]),
# c(mx, bbox[3], my, bbox[4])
# )
#
# var_name0 = rbind(
# c('soc', '_0-5cm_mean'),
# c('soc', '_5-15cm_mean'),
# c('soc', '_15-30cm_mean'),
# c('phh2o', '_0-5cm_mean'),
# c('phh2o', '_5-15cm_mean'),
# c('phh2o', '_15-30cm_mean')
# )
#
# ras_tot = list()
# for(vp in nrow(var_name0)){
# var_name3 = paste0(var_name0[vp,], collapse = '')
# cat('\n', var_name3, '\n\t')
#
# part1_url = paste0('https://maps.isric.org/mapserv?map=/map/',
# var_name0[vp,1],
# '.map&SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=',
# var_name3,
# '&FORMAT=image/tiff&SUBSET=')
# part3_url = '&SUBSETTINGCRS=http://www.opengis.net/def/crs/EPSG/0/4326&OUTPUTCRS=http://www.opengis.net/def/crs/EPSG/0/4326'
# ras = list()
# for(i in 1:4){
# cat(i, '')
# part2_url = f_part2_url(ras_coor[i,])
#
# full_url = paste0(part1_url, part2_url, part3_url)
# path_ras = paste0('datos/',var_name3,'_',i,'.tif')
# r = gdal_translate(full_url, path_ras)
#
# ras[[i]] = raster(path_ras)
# }
# ras_tot[[var_name3]] = merge(ras[[1]],ras[[2]],ras[[3]],ras[[4]])
# writeRaster(ras_tot[[var_name3]],
# paste0('datos/merge/',var_name3,'_tot.tif'),
# overwrite=TRUE)
# }
# names(ras_tot)
Unir los raster sumando de cada profundidad
- Cargar las 3 profundidades
- Para el COS total se suman las capas, ya que estas estan en dg (decigramos) se divide entre 10000 para tenerlos en \(COS~kg/kg_{suelo}\)
- Para el pH se calculo el promedio en toda la profundidad, además se dividio entre 10 ya que los datos provienen multiplicados por este valor
# base_path = '../Clase 7 - Suelo - 161221/datos/merge/'
# vari = c('soc', 'phh2o')
# depths = c('0-5', '5-15', '15-30')
#
# soc_0_5 = raster(paste0(base_path, 'soc_', depths[1], 'cm_mean_tot.tif'))
# soc_5_15 = raster(paste0(base_path, 'soc_', depths[2], 'cm_mean_tot.tif'))
# soc_15_30 = raster(paste0(base_path, 'soc_', depths[3], 'cm_mean_tot.tif'))
# soc_0_30 = (soc_0_5 + soc_5_15 + soc_15_30)/(10*1000)
# writeRaster(soc_0_30, paste0('merge_0_30/soc_0_30cm.tif'))
#
# phh2o_0_5 = raster(paste0(base_path, 'phh2o_', depths[1], 'cm_mean_tot.tif'))
# phh2o_5_15 = raster(paste0(base_path, 'phh2o_', depths[2], 'cm_mean_tot.tif'))
# phh2o_15_30 = raster(paste0(base_path, 'phh2o_', depths[3], 'cm_mean_tot.tif'))
# phh2o_0_30 = (phh2o_0_5 + phh2o_5_15 + phh2o_15_30)/(10*3)
# writeRaster(phh2o_0_30, paste0('merge_0_30/phh2o_0_30cm.tif'))
Capas finales
soc = raster('merge_0_30/soc_0_30cm.tif')
ph = raster('merge_0_30/phh2o_0_30cm.tif')
SOC
soc_aru = mask(soc, shp_aru)
soc_aru_pto = rasterToPoints(soc_aru, spatial = T)
soc_aru_magna = st_transform(st_as_sf(soc_aru_pto), crs=3116)
soc_aru_magna2 = as(soc_aru_magna, 'Spatial')
tm_shape(soc_aru)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Aumento de escala
soc_agg = aggregate(soc_aru, 30)
soc_agg = mask(soc_agg, shp_aru)
# writeRaster(soc_agg, 'soc_agg.tif')
soc_agg_pto = rasterToPoints(soc_agg, spatial = T)
soc_agg_df = as.data.frame(soc_agg_pto)
soc_magna = st_transform(st_as_sf(soc_agg_pto), crs=3116)
soc_magna2 = as(soc_magna, 'Spatial')
tm_shape(soc_agg)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Muestreo
set.seed(123)
(corte_tt = floor(0.0005*nrow(soc_aru_magna2)))
## [1] 179
index30 = sample(nrow(soc_aru_magna2), corte_tt)
soc_test = soc_aru_magna2[index30, ]
plot(soc_magna2)
points(soc_test, pch = 16)
PolÃgonos de Thiessen
vor = voronoi_polygon(soc_agg_df)
crs(vor) = crs(soc_agg_pto)
vor_df = fortify_voronoi(vor)
ggplot(vor_df) +
geom_polygon(aes(x=x, y=y ,fill=point.x, group=group),
size=0, color = gray(0.5))+
scale_fill_distiller(palette = 'Oranges', direction = 1)+
geom_sf(data=shp_aru, fill='transparent', col = 'black')+
theme_bw()
IDW
powers = seq(3, 3.1, .001)
rmse_idw = NULL
for(p in powers){
idw_p = gstat::idw(soc_0_30cm ~ 1, soc_magna2, newdata=soc_test, idp=p)
rmse_idw = c(rmse_idw, Metrics::rmse(soc_test$soc_0_30cm, idw_p$var1.pred))
}
(bp = powers[which.min(rmse_idw)])
## [1] 3.063
ggplot()+
geom_line(aes(powers, rmse_idw), size = 1)+
geom_point(aes(bp, min(rmse_idw)), size = 2, col = 'blue')+
theme_bw()
grd = as.data.frame(spsample(soc_magna2, "regular", n=10000))
names(grd) = c("X", "Y")
coordinates(grd) = c("X", "Y")
gridded(grd) = TRUE
fullgrid(grd) = TRUE
crs(grd) = crs(soc_magna2)
cord_test = coordinates(soc_test)
cord_test = as.data.frame(cord_test)
colnames(cord_test) = c('X', 'Y')
soc_idw = gstat::idw(soc_0_30cm ~ 1, soc_magna2, newdata=grd, idp=bp)
## [inverse distance weighted interpolation]
soc_idw_r = raster(soc_idw)
soc_idw_r_aru = raster::mask(soc_idw_r, shp_aru_magna)
# rmse
idw_p = gstat::idw(soc_0_30cm ~ 1, soc_magna2, newdata=soc_test, idp=p)
## [inverse distance weighted interpolation]
rmse_idw = Metrics::rmse(soc_test$soc_0_30cm, idw_p$var1.pred)
ggplot()+
geom_raster(data = as.data.frame(soc_idw_r_aru, xy=T),
aes(x, y, fill = var1.pred)) +
scale_fill_gradientn(colors = brewer.pal(n=8, name = 'Oranges'),
na.value="transparent",
name = 'SOC fitted IDW')+
geom_sf(data = shp_aru_magna, fill=NA, size=1)+
geom_sf_text(data = shp_aru_magna, aes(label = MPIO_CNMBR), size = 3)+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'),
text = element_text(size = 15))
tm_shape(soc_idw_r_aru)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Ajuste polinomial de 1er orden
P = soc_magna2
f.1 = as.formula(soc_0_30cm ~ X + Y)
P$X = coordinates(P)[,1]
P$Y = coordinates(P)[,2]
lm.1 = lm(f.1, data=P)
# rmse
pred_test = predict(lm.1, cord_test)
rmse_p1 = Metrics::rmse(soc_test$soc_0_30cm, pred_test)
dat.1st = SpatialGridDataFrame(
grd,
data.frame(var1.pred = predict(lm.1,newdata=grd))
)
r = raster(dat.1st)
r.m = raster::mask(r, shp_aru_magna)
tm_shape(r.m)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Ajuste polinomial de 2do orden
f.2 = as.formula(soc_0_30cm ~ X + Y + I(X*X)+I(Y*Y) + I(X*Y))
P$X = coordinates(P)[,1]
P$Y = coordinates(P)[,2]
lm.2 = lm(f.2, data=P)
# rmse
pred_test = predict(lm.2, cord_test)
rmse_p2 = Metrics::rmse(soc_test$soc_0_30cm, pred_test)
dat.2nd = SpatialGridDataFrame(
grd,
data.frame(var1.pred = predict(lm.2, newdata=grd))
)
r = raster(dat.2nd)
r.m = raster::mask(r, shp_aru_magna)
tm_shape(r.m)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Kriging
f.1 = as.formula(soc_0_30cm ~ X + Y)
# var.smpl = variogram(f.1, P, cloud = FALSE, cutoff=200000, width=10000)
var.smpl = variogram(f.1, P)
f = function(x) {
attr(m.fit <<- fit.variogram(
var.smpl,
vgm("Gau", nugget=x,kappa=NULL)), "SSErr")
}
optimize(f, seq(0.1, 5, 0.1))
## $minimum
## [1] 4.999924
##
## $objective
## [1] 6.133919e-13
m.fit
## model psill range
## 1 Nug 2.643084e-05 0.00
## 2 Gau 1.674362e-03 21846.98
# dat.fit = fit.variogram(
# var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
# vgm(psill=5, model="Gau", range=120000, nugget=-0.5)
# )
dat.fit = fit.variogram(var.smpl, vgm(model = 'Gau',
nugget = 18620.27,
psill = 148893.42,
range = 22194.62))
dat.fit
## model psill range
## 1 Nug 0.0000397995 0.00
## 2 Gau 0.0016587222 22912.54
plot(var.smpl, dat.fit)
grd_krg = grd
dat.krg = krige(f.1, P, grd_krg, dat.fit)
## [using universal kriging]
r = raster(dat.krg)
r.m = raster::mask(r, shp_aru_magna)
tm_shape(r.m)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
AKIMA
a = interp(P$X,
P$Y,
P$soc_0_30cm,
nx = 1000, ny = 1000,
linear = F)
r = raster(a)
r.m = raster::mask(r, shp_aru_magna)
crs(r.m) = crs(shp_aru_magna)
tm_shape(r.m)+
tm_raster(n = 5)+
tm_shape(shp_aru_magna)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
RMSE
rmse_idw
## [1] 0.01750904
rmse_p1
## [1] 0.04246401
rmse_p2
## [1] 0.03299804
pH
ph_aru = mask(ph, shp_aru)
ph_aru_pto = rasterToPoints(ph_aru, spatial = T)
ph_aru_magna = st_transform(st_as_sf(ph_aru_pto), crs=3116)
ph_aru_magna2 = as(ph_aru_magna, 'Spatial')
tm_shape(ph_aru)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Aumento de escala
ph_agg = aggregate(ph_aru, 30)
ph_agg = mask(ph_agg, shp_aru)
# writeRaster(ph_agg, 'ph_agg.tif')
ph_agg_pto = rasterToPoints(ph_agg, spatial = T)
ph_agg_df = as.data.frame(ph_agg_pto)
ph_magna = st_transform(st_as_sf(ph_agg_pto), crs=3116)
ph_magna2 = as(ph_magna, 'Spatial')
tm_shape(ph_agg)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Muestreo
set.seed(123)
(corte_tt = floor(0.0005*nrow(ph_aru_magna2)))
## [1] 179
index30 = sample(nrow(ph_aru_magna2), corte_tt)
ph_test = ph_aru_magna2[index30, ]
plot(ph_magna2)
points(ph_test, pch = 16)
PolÃgonos de Thiessen
vor = voronoi_polygon(ph_agg_df)
crs(vor) = crs(ph_agg_pto)
vor_df = fortify_voronoi(vor)
ggplot(vor_df) +
geom_polygon(aes(x=x, y=y ,fill=point.x, group=group),
size=0, color = gray(0.5))+
scale_fill_distiller(palette = 'Oranges', direction = 1)+
geom_sf(data=shp_aru, fill='transparent', col = 'black')+
theme_bw()
IDW
powers = seq(3.6, 3.8, .01)
rmse_idw = NULL
for(p in powers){
idw_p = gstat::idw(phh2o_0_30cm ~ 1, ph_magna2, newdata=ph_test, idp=p)
rmse_idw = c(rmse_idw, Metrics::rmse(ph_test$phh2o_0_30cm, idw_p$var1.pred))
}
(bp = powers[which.min(rmse_idw)])
## [1] 3.67
ggplot()+
geom_line(aes(powers, rmse_idw), size = 1)+
geom_point(aes(bp, min(rmse_idw)), size = 2, col = 'blue')+
theme_bw()
grd = as.data.frame(spsample(ph_magna2, "regular", n=10000))
names(grd) = c("X", "Y")
coordinates(grd) = c("X", "Y")
gridded(grd) = TRUE
fullgrid(grd) = TRUE
crs(grd) = crs(ph_magna2)
ph_idw = gstat::idw(phh2o_0_30cm ~ 1, ph_magna2, newdata=grd, idp=bp)
## [inverse distance weighted interpolation]
ph_idw_r = raster(ph_idw)
ph_idw_r_aru = raster::mask(ph_idw_r, shp_aru_magna)
ggplot()+
geom_raster(data = as.data.frame(ph_idw_r_aru, xy=T),
aes(x, y, fill = var1.pred)) +
scale_fill_gradientn(colors = brewer.pal(n=8, name = 'Oranges'),
na.value="transparent",
name = 'ph fitted IDW')+
geom_sf(data = shp_aru_magna, fill=NA, size=1)+
geom_sf_text(data = shp_aru_magna, aes(label = MPIO_CNMBR), size = 3)+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'),
text = element_text(size = 15))
tm_shape(ph_idw_r_aru)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Ajuste polinomial de 1er orden
P = ph_magna2
f.1 = as.formula(phh2o_0_30cm ~ X + Y)
P$X = coordinates(P)[,1]
P$Y = coordinates(P)[,2]
lm.1 = lm(f.1, data=P)
dat.1st = SpatialGridDataFrame(
grd,
data.frame(var1.pred = predict(lm.1,newdata=grd))
)
r = raster(dat.1st)
r.m = raster::mask(r, shp_aru_magna)
tm_shape(r.m)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Ajuste polinomial de 2do orden
f.2 = as.formula(phh2o_0_30cm ~ X + Y + I(X*X)+I(Y*Y) + I(X*Y))
P$X = coordinates(P)[,1]
P$Y = coordinates(P)[,2]
lm.2 = lm(f.2, data=P)
dat.2nd = SpatialGridDataFrame(
grd,
data.frame(var1.pred = predict(lm.2, newdata=grd))
)
r = raster(dat.2nd)
r.m = raster::mask(r, shp_aru_magna)
tm_shape(r.m)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
Kriging
f.1 = as.formula(phh2o_0_30cm ~ X + Y)
# var.smpl = variogram(f.1, P, cloud = FALSE, cutoff=200000, width=10000)
var.smpl = variogram(f.1, P)
f = function(x) {
attr(m.fit <<- fit.variogram(
var.smpl,
vgm("Gau", nugget=x,kappa=NULL)), "SSErr")
}
optimize(f, seq(0.1, 5, 0.1))
## $minimum
## [1] 4.999924
##
## $objective
## [1] 4.28698e-11
m.fit
## model psill range
## 1 Nug 0.006820373 0.0
## 2 Gau 0.012006054 17542.3
dat.fit = fit.variogram(var.smpl, vgm(model = 'Gau',
nugget = 33.01304,
psill = 68.25435,
range = 24149.93))
dat.fit
## model psill range
## 1 Nug 0.004663404 0
## 2 Gau 0.014149658 11577
plot(var.smpl, dat.fit)
grd_krg = grd
dat.krg = krige(f.1, P, grd_krg, dat.fit)
## [using universal kriging]
r = raster(dat.krg)
r.m = raster::mask(r, shp_aru_magna)
tm_shape(r.m)+
tm_raster(n=5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
AKIMA
a = interp(P$X,
P$Y,
P$phh2o_0_30cm,
nx = 1000, ny = 1000,
linear = F)
r = raster(a)
r.m = raster::mask(r, shp_aru_magna)
crs(r.m) = crs(shp_aru_magna)
tm_shape(r.m)+
tm_raster(n = 5)+
tm_shape(shp_aru)+
tm_polygons(alpha = 0, lwd = 3, border.col = 'black')+
tm_legend(outside = TRUE)+
tm_grid(n.x = 5, n.y = 5, col = gray(0.7))
