El IDW es un tipo de mĆ©todo determinista para la interpolaciĆ³n espacial con un conjunto de puntos conocidos y dispersos (Chen et al., 2014). Su idea general se basa en la suposiciĆ³n de que el valor del atributo de un punto no muestreado es la media ponderada de los valores conocidos dentro de la vecindad, y los pesos estĆ”n relacionados inversamente con las distancias entre la ubicaciĆ³n de la predicciĆ³n y las ubicaciones muestreadas (Lu & Wong, 2008).
rm(list=ls())library(sf)Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUElibrary(rgdal)Loading required package: sp
Please note that rgdal will be retired by the end of 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
rgdal: version: 1.5-32, (SVN revision 1176)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
Path to GDAL shared files: C:/Users/ynata/AppData/Local/R/win-library/4.2/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
Path to PROJ shared files: C:/Users/ynata/AppData/Local/R/win-library/4.2/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.5-0
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(rasterVis)Loading required package: latticelibrary(dplyr)
Attaching package: ādplyrā
The following objects are masked from āpackage:statsā:
filter, lag
The following objects are masked from āpackage:baseā:
intersect, setdiff, setequal, unionlibrary(raster)
Attaching package: ārasterā
The following object is masked from āpackage:dplyrā:
selectlibrary(gstat)
#library(geoR)
library(berryFunctions)
Attaching package: āberryFunctionsā
The following object is masked from āpackage:rasterā:
distance
The following object is masked from āpackage:dplyrā:
betweenmun.tmp = st_read("C:\\Users\\ynata\\OneDrive\\Documentos\\GGB2022\\Datos\\DPTOCASANARE.shp")Reading layer `DPTOCASANARE' from data source
`C:\Users\ynata\OneDrive\Documentos\GGB2022\Datos\DPTOCASANARE.shp'
using driver `ESRI Shapefile'
Simple feature collection with 19 features and 11 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
Geodetic CRS: WGS 84mun.tmp %>% dplyr::select( MPIO_CCNCT, MPIO_CNMBR) -> casanare
rename(casanare, MPIO_CCDGO = MPIO_CCNCT)Simple feature collection with 19 features and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
Geodetic CRS: WGS 84
First 10 features:
MPIO_CCDGO MPIO_CNMBR geometry
1 85001 YOPAL POLYGON ((-72.39513 5.56853...
2 85010 AGUAZUL POLYGON ((-72.56545 5.36972...
3 85015 CHĆMEZA POLYGON ((-72.81017 5.36659...
4 85136 LA SALINA POLYGON ((-72.33885 6.34470...
5 85139 MANĆ POLYGON ((-72.34155 5.06495...
6 85162 MONTERREY POLYGON ((-72.89989 5.03465...
7 85225 NUNCHĆA POLYGON ((-72.19558 5.71924...
8 85230 OROCUĆ POLYGON ((-71.50965 5.21730...
9 85263 PORE POLYGON ((-72.04587 5.83819...
10 85279 RECETOR POLYGON ((-72.80501 5.38092...#file.choose()precip <- raster("C:\\Users\\ynata\\OneDrive\\Documentos\\GGB2022\\R\\precip.tif")precipclass : RasterLayer
dimensions : 51, 78, 3978 (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -73.08333, -69.83333, 4.25, 6.375 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : precip.tif
names : precip
values : 0.2, 161.1 (min, max)## change attribute name
names(precip) <- "X2019.02"library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "#fcc000","yellow", "cyan", "blue", "#3240cd"), values(precip$X2019.02),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(precip$X2019.02 , colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(precip$X2019.02),
title = "Rainfall-Feb.2019 [mm]")(precip.mask <- mask(x = precip, mask = casanare))class : RasterLayer
dimensions : 51, 78, 3978 (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -73.08333, -69.83333, 4.25, 6.375 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : X2019.02
values : 0.2, 122.2 (min, max)leaflet() %>% addTiles() %>%
addRasterImage(precip.mask$X2019.02 , colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(precip.mask$X2019.02),
title = "Rainfall in Casanare-Feb.2019 [mm]")(precip.points <- rasterToPoints(precip.mask, spatial = TRUE))class : SpatialPointsDataFrame
features : 2086
extent : -73.0625, -69.85417, 4.3125, 6.270833 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 1
names : X2019.02
min values : 0.200000002980232
max values : 122.199996948242 names(precip.points) <- "rain"# crs=32618
(p.sf.magna <- st_transform(st_as_sf(precip.points), crs=32618))Simple feature collection with 2086 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 714880.9 ymin: 477045.6 xmax: 1070220 ymax: 695592.3
Projected CRS: WGS 84 / UTM zone 18N
First 10 features:
rain geometry
1 26.8 POINT (788153.7 693858.9)
2 0.5 POINT (912783 694612.1)
3 0.5 POINT (917401.6 694645.1)
4 0.5 POINT (922020.5 694678.6)
5 1.1 POINT (1005199 695343.8)
6 1.0 POINT (1009823 695384.3)
7 0.9 POINT (1014447 695425.2)
8 1.3 POINT (1019071 695466.4)
9 1.1 POINT (1023695 695508)
10 1.0 POINT (1028319 695550)(casanare.sf.magna <- st_transform(casanare, crs=32618))Simple feature collection with 19 features and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 713187.7 ymin: 474289.2 xmax: 1072227 ymax: 702220.4
Projected CRS: WGS 84 / UTM zone 18N
First 10 features:
MPIO_CCNCT MPIO_CNMBR geometry
1 85001 YOPAL POLYGON ((788595.8 616146.7...
2 85010 AGUAZUL POLYGON ((769803.5 594069.8...
3 85015 CHĆMEZA POLYGON ((742671.1 593621.7...
4 85136 LA SALINA POLYGON ((794421.3 702065.3...
5 85139 MANĆ POLYGON ((794776.7 560448.2...
6 85162 MONTERREY POLYGON ((732846.2 556869.8...
7 85225 NUNCHĆA POLYGON ((810641.1 632927, ...
8 85230 OROCUĆ POLYGON ((887027.6 577757.7...
9 85263 PORE POLYGON ((827170.2 646176.4...
10 85279 RECETOR POLYGON ((743237.1 595208.9...precip2 <- as(p.sf.magna, 'Spatial')casanare2 <- as(casanare.sf.magna, 'Spatial')precip2@bbox min max
coords.x1 714880.9 1070220.1
coords.x2 477045.6 695592.3casanare2@bbox min max
x 713187.7 1072227.0
y 474289.2 702220.4precip2@bbox <- casanare2@bboxtrain_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,]ptos_trainclass : SpatialPointsDataFrame
features : 1460
extent : 719449.4, 1070220, 477045.6, 695592.3 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
variables : 1
names : rain
min values : 0.200000002980232
max values : 122.199996948242 library(gstat)
# Crear una cuadrĆcula vacĆa donde n es el nĆŗmero total de celdas
grd <- as.data.frame(spsample(precip2, "regular", n=319000))
# Tienes que averiguar cuƔl es el tamaƱo esperado del grd de salida
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE # Crear un objeto SpatialPixel
fullgrid(grd) <- TRUE # Crear el objeto SpatialGrid# Definir el CRS mediante un cĆ³digo EPSG
x <- CRS(SRS_string='EPSG:32618')
# Mostrar el CRS almacenado usando comment()
#cat(comment(x), "\n")
# Almacenar el wkt en una variable
wkt <- comment(x)
# Usa esto para asignar el CRS de otro objeto sp
(y <- CRS(SRS_string = wkt))Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["WGS 84 / UTM zone 18N",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 18N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Between 78Ā°W and 72Ā°W, northern hemisphere between equator and 84Ā°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela."],
BBOX[0,-78,84,-72]],
ID["EPSG",32618]] # AƱadir la informaciĆ³n de proyecciĆ³n de P a la rejilla vacĆa
sp::proj4string(grd) <- y# Interpolar las celdas de la cuadrĆcula utilizando un valor de potencia de 2 (idp=2,0)
P.idw <- gstat::idw(rain ~ 1, ptos_train, newdata=grd, idp=2.0)[inverse distance weighted interpolation]# Convertir a objeto rasterizado y luego recortar a AOI
r <- raster(P.idw)(r.m <- raster::mask(r, casanare2))class : RasterLayer
dimensions : 450, 708, 318600 (nrow, ncol, ncell)
resolution : 506.4978, 506.4978 (x, y)
extent : 713396.2, 1071997, 474137.9, 702061.9 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
source : memory
names : var1.pred
values : 0.3206528, 120.5286 (min, max)leaflet() %>% addTiles() %>%
addRasterImage(r.m$var1.pred, colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(r.m$var1.pred),
title = "IDW rainfall in February 2019 [mm]")