Códigos para la interpolación

Polinomio de Primer orden

precip.points <- geojsonio::geojson_read("./chirps/ppoints.geojson", what="sp")

(Vichada <- shapefile("./ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 4 
extent      : -71.07793, -67.4098, 2.737109, 6.324317  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO,    MPIO_CNMBR,                           MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,     Shape_Area 
min values  :         99,      99001,      CUMARIBO,       Decreto 1594 de Ago 5  de 1974, 3898.56891769,      2017,    VICHADA, 3.29670807195, 0.316732435778 
max values  :         99,      99773, SANTA ROSALÍA, Ordenanza 66 de Noviembre 22 de 1996, 65599.7022767,      2017,    VICHADA,  18.794382661,   5.3085802966 

Vichada_sf <-  sf::st_as_sf(Vichada)
(border_sf <- Vichada_sf %>%summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -71.07793 ymin: 2.737109 xmax: -67.4098 ymax: 6.324317
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
      area                       geometry
1 100065.9 POLYGON ((-67.814 5.342584,...

 (border <- as(border_sf, 'Spatial')) 

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -71.07793, -67.4098, 2.737109, 6.324317  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       :            area 
value       : 100065.86486549 

(Vichada.sf <- st_as_sf(Vichada) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))

Simple feature collection with 4 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -71.07793 ymin: 2.737109 xmax: -67.4098 ymax: 6.324317
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
           MUNIC CODIGO
1  SANTA ROSALÍA  99624
2 PUERTO CARREÑO  99001
3   LA PRIMAVERA  99524
4       CUMARIBO  99773
                        geometry
1 POLYGON ((-70.65378 5.37297...
2 POLYGON ((-67.80972 6.32431...
3 POLYGON ((-69.03359 6.21869...
4 POLYGON ((-68.47074 5.55046...

p.sf <- st_as_sf(precip.points)
(precip.sf = st_intersection(Vichada.sf, p.sf))

although coordinates are longitude/latitude, st_intersection assumes that they are planar
attribute variables are assumed to be spatially constant throughout all geometries

Simple feature collection with 3248 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -71.075 ymin: 2.774999 xmax: -67.475 ymax: 6.274999
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
First 10 features:
             MUNIC CODIGO    lluvia
2   PUERTO CARREÑO  99001 0.8161542
2.1 PUERTO CARREÑO  99001 0.8258479
2.2 PUERTO CARREÑO  99001 0.8074083
2.3 PUERTO CARREÑO  99001 0.8239177
2.4 PUERTO CARREÑO  99001 0.7755996
2.5 PUERTO CARREÑO  99001 0.7823746
2.6 PUERTO CARREÑO  99001 3.5164738
2.7 PUERTO CARREÑO  99001 3.2473364
2.8 PUERTO CARREÑO  99001 0.7987307
2.9 PUERTO CARREÑO  99001 0.7961798
                    geometry
2   POINT (-67.875 6.274999)
2.1 POINT (-67.825 6.274999)
2.2 POINT (-67.775 6.274999)
2.3 POINT (-67.725 6.274999)
2.4 POINT (-67.675 6.274999)
2.5 POINT (-67.625 6.274999)
2.6 POINT (-68.125 6.224999)
2.7 POINT (-68.075 6.224999)
2.8 POINT (-67.925 6.224999)
2.9 POINT (-67.875 6.224999)

p.sf.magna <- st_transform(precip.sf, crs=3116)
Vichada.sf.magna <- st_transform(Vichada.sf, crs=3116)
(precip2 <- as(p.sf.magna, 'Spatial'))

class       : SpatialPointsDataFrame 
features    : 3248 
extent      : 1333188, 1732620, 799126.6, 1190067  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +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 
variables   : 3
names       :         MUNIC, CODIGO,            lluvia 
min values  :      CUMARIBO,  99001, 0.541795015335083 
max values  : SANTA ROSALÍA,  99773,  85.8214721679688 

shapefile(precip2, filename='./chirps/precip2.shp', overwrite=TRUE)

precip2$precipitacion <- round(precip2$lluvia, 1)
(Vichada2 <- as(Vichada.sf.magna, 'Spatial'))

class       : SpatialPolygonsDataFrame 
features    : 4 
extent      : 1332850, 1739856, 794950.2, 1195302  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +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 
variables   : 2
names       :         MUNIC, CODIGO 
min values  :      CUMARIBO,  99001 
max values  : SANTA ROSALÍA,  99773 

shapefile(Vichada2, filename='./chirps/Vichada2.shp', overwrite=TRUE)
precip2@bbox <- Vichada2@bbox

train_index <- sample(1:nrow(precip2), 0.7 * nrow(precip2))
test_index <- setdiff(1:nrow(precip2), train_index)
ptos_train <- precip2[train_index, ]
ptos_test  <- precip2[test_index,]

P <- ptos_train

f.1 <- as.formula(precipitacion ~ X + Y) 

P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]

lm.1 <- lm( f.1, data=P)

grd              <- as.data.frame(spsample(ptos_train, "regular", n=500000))

names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE 
fullgrid(grd)    <- TRUE  
proj4string(grd) <- proj4string(ptos_train)

dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 

r   <- raster(dat.1st)
r.m1 <- raster::mask(r, Vichada2)

tm_shape(r.m1) + 
  tm_raster(n=10, palette="RdBu", midpoint = 43, 
            title="Interpolación con polinomio de 1er orden \n para precipitacion \n(mm)") +
  tm_shape(P) + tm_dots(size=0.02) +
  tm_legend(legend.outside=TRUE)

2.KRIGING

f.1 <- as.formula(precipitacion ~ X + Y) 

2.1 Variograma Experimental

var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=200000, width=10000)

plot(var.smpl)

2.2Variograma de Ajuste a los Datos

Ensayo <-  autofitVariogram(f.1, P)
summary(Ensayo)

Experimental variogram:

Fitted variogram model:

Sums of squares betw. var. model and sample var.[1] 0.003434569

plot(Ensayo, pch=16, col="black")

dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE, vgm(psill=68.711235, model="Gau", range=46602.27   , nugget=4))

plot(var.smpl, dat.fit, xlim=c(0,150000))

2.3 Superficie de Interpolación

f.1 <- as.formula(precipitacion ~ X + Y) 

dat.krg <- krige( f.1, P, grd, dat.fit)

CRS object has comment, which is lost in output

[using universal kriging]

r1 <- raster(dat.krg)
r.k <- raster::mask(r1, Vichada2)
r.k

class      : RasterLayer 
dimensions : 695, 720, 500400  (nrow, ncol, ncell)
resolution : 554.8374, 554.8374  (x, y)
extent     : 1333051, 1732534, 804470.7, 1190083  (xmin, xmax, ymin, ymax)
crs        : +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 
source     : memory
names      : var1.pred 
values     : 0.1595495, 80.58832  (min, max)

tm_shape(r.k) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, midpoint = NA,
            title="Universal Kriging\npara precipitacion \n(mm)") +
  tm_shape(P) + tm_dots(size=0.02, col="white") +
  tm_legend(legend.outside=TRUE)

The argument auto.palette.mapping is deprecated. Please use midpoint for numeric data and stretch.palette for categorical data to control the palette mapping.

r2   <- raster(dat.krg, layer="var1.var")
r.j <- raster::mask(r2, Vichada2)

tm_shape(r.j) + 
  tm_raster(n=7, palette ="Reds",
            title="Kriging Interpolation\nVariance map \n(in squared mm)") +tm_shape(P) + tm_dots(size=0.02) +
  tm_legend(legend.outside=TRUE)

r3   <- sqrt(raster(dat.krg, layer="var1.var")) * 1.96
r.l <- raster::mask(r3, Vichada2)

tm_shape(r.l) + 
  tm_raster(n=7, palette ="Reds",
            title="Kriging Interpolation\n95% CI map \n(in mm)") +tm_shape(P) + tm_dots(size=0.02, col="white") +
  tm_legend(legend.outside=TRUE)

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8    
 [5] LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C             
 [9] LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] tmap_2.3-1   gstat_2.0-6  raster_3.1-5 sp_1.4-2    

loaded via a namespace (and not attached):
 [1] zoo_1.8-8          xfun_0.14          sf_0.9-3           lattice_0.20-38    V8_3.1.0          
 [6] htmltools_0.5.0    viridisLite_0.3.0  spacetime_1.2-3    yaml_2.2.1         base64enc_0.1-3   
[11] XML_3.99-0.3       rlang_0.4.6        e1071_1.7-3        foreign_0.8-75     httpcode_0.3.0    
[16] DBI_1.1.0          RColorBrewer_1.1-2 stringr_1.4.0      geojsonio_0.9.2    rgeos_0.5-3       
[21] htmlwidgets_1.5.1  leafsync_0.1.0     codetools_0.2-16   evaluate_0.14      knitr_1.28        
[26] maptools_1.0-1     crosstalk_1.1.0.1  curl_4.3           class_7.3-15       xts_0.12-0        
[31] Rcpp_1.0.4.6       KernSmooth_2.23-16 classInt_0.4-3     lwgeom_0.2-1       leaflet_2.0.3     
[36] jsonlite_1.6.1     FNN_1.1.3          packrat_0.5.0      digest_0.6.25      stringi_1.4.6     
[41] tmaptools_3.0      grid_3.6.3         jqr_1.1.0          rgdal_1.4-8        tools_3.6.3       
[46] magrittr_1.5       lazyeval_0.2.2     geojson_0.3.2      dichromat_2.0-0    crul_0.9.0        
[51] rsconnect_0.8.16   rmarkdown_2.1      R6_2.4.1           units_0.6-6        intervals_0.15.2  
[56] compiler_3.6.3