Interpolación - datos de precipitación de La Guajira

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(knitr)
#include_graphics('./chirps/chirps_pentad.png')
library(knitr)
#include_graphics('./chirps/chirps_pentad_2020.png')
library(rgdal)

Loading required package: sp
rgdal: version: 1.5-10, (SVN revision 1006)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
Path to GDAL shared files: C:/Users/familiar/Documents/R/R-3.6.3/library/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: C:/Users/familiar/Documents/R/R-3.6.3/library/rgdal/proj
Linking to sp version:1.4-2
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

## Loading required package: sp
## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-2
library(raster)
library(sf)

Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.2     v purrr   0.3.4
v tibble  3.0.1     v dplyr   1.0.0
v tidyr   1.1.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x tidyr::extract() masks raster::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x dplyr::select()  masks raster::select()

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
# ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.3     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
# ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()
library(tmap)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(gstat)
library(sp)

## Read the CHIRPS dataset
precip <- raster("C:/Users/familiar/Documents/Geo R final/chirps-v2.0.2020.05.6.tif")

(aoi <- shapefile("C:/Users/familiar/Downloads/gb2/La guajira/44_LA_GUAJIRA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 15 
extent      : -73.66494, -71.11296, 10.39676, 12.45944  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                               MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,     Shape_Leng,      Shape_Area 
min values  :         44,      44001,    ALBANIA,                                     1888,  179.59085022,      2017, LA GUAJIRA, 0.710049632727, 0.0148255649224 
max values  :         44,      44874, VILLANUEVA, Ordenanza015 del 29 de Noviembre de 1989, 7886.05230923,      2017, LA GUAJIRA,  6.11271226506,  0.653549310217 

precip.crop <- raster::crop(precip, extent(aoi))
precip.mask <- mask(x = precip.crop, mask = aoi)
precip.mask

class      : RasterLayer 
dimensions : 41, 51, 2091  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -73.65, -71.1, 10.4, 12.45  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : chirps.v2.0.2020.05.6 
values     : 0, 177.7225  (min, max)

plot(precip.mask, main= "CHIRPS Lluvias en La Guajira desde el  01.06  al  05.06 en 2020 [mm]")
plot(aoi, add=TRUE)

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("orange", "yellow", "green3","cyan", "blue", "darkblue"), values(precip.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(precip.mask, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(precip.mask),
    title = "CHIRPS Lluvias en La Guajira desde el  01.06 al 05.06 en 2020 [mm]")

Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsSome values were outside the color scale and will be treated as NA

precip.points <- rasterToPoints(precip.mask, spatial = TRUE)
precip.points

class       : SpatialPointsDataFrame 
features    : 684 
extent      : -73.625, -71.125, 10.425, 12.425  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       : chirps.v2.0.2020.05.6 
min values  :                     0 
max values  :      177.722457885742 

names(precip.points) <- "rain"
precip.points

class       : SpatialPointsDataFrame 
features    : 684 
extent      : -73.625, -71.125, 10.425, 12.425  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :             rain 
min values  :                0 
max values  : 177.722457885742 

str(precip.points)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 684 obs. of  1 variable:
  .. ..$ rain: num [1:684] 2.92 4.17e-01 6.77e-15 0.00 1.54 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:684, 1:2] -71.7 -71.6 -71.6 -71.5 -71.7 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -73.6 10.4 -71.1 12.4
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"

plot(precip.mask, main= "CHIRPS Lluvias en La Guajira desde el  01.06 al 05.06 en 2020 [mm]")
plot(aoi, add=TRUE)

points(precip.points$x, precip.points$y, col = "red", cex = .6)

library(geojsonio)
geojsonio::geojson_write(precip.points, file = "./chirps/ppoints.geojson")

Success! File is at ./chirps/ppoints.geojson

<geojson-file>
  Path:       ./chirps/ppoints.geojson
  From class: SpatialPointsDataFrame

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

class       : SpatialPointsDataFrame 
features    : 684 
extent      : -73.625, -71.125, 10.425, 12.425  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :             rain 
min values  :                0 
max values  : 177.722457885742 

(aoi <- shapefile("C:/Users/familiar/Downloads/gb2/La guajira/44_LA_GUAJIRA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 15 
extent      : -73.66494, -71.11296, 10.39676, 12.45944  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                               MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,     Shape_Leng,      Shape_Area 
min values  :         44,      44001,    ALBANIA,                                     1888,  179.59085022,      2017, LA GUAJIRA, 0.710049632727, 0.0148255649224 
max values  :         44,      44874, VILLANUEVA, Ordenanza015 del 29 de Noviembre de 1989, 7886.05230923,      2017, LA GUAJIRA,  6.11271226506,  0.653549310217 

guajira_sf <-  sf::st_as_sf(aoi)
(border_sf <-
  guajira_sf %>%
  summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -73.66494 ymin: 10.39676 xmax: -71.11296 ymax: 12.45944
geographic CRS: WGS 84
      area                       geometry
1 20621.92 POLYGON ((-72.88683 10.4532...

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

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -73.66494, -71.11296, 10.39676, 12.45944  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :           area 
value       : 20621.92459446 

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

Simple feature collection with 15 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.66494 ymin: 10.39676 xmax: -71.11296 ymax: 12.45944
geographic CRS: WGS 84
First 10 features:
                MUNIC CODIGO                       geometry
1           BARRANCAS  44078 MULTIPOLYGON (((-72.59885 1...
2        DISTRACCIÓN  44098 MULTIPOLYGON (((-72.89246 1...
3           EL MOLINO  44110 MULTIPOLYGON (((-72.96901 1...
4             FONSECA  44279 MULTIPOLYGON (((-72.91056 1...
5           HATONUEVO  44378 MULTIPOLYGON (((-72.6568 11...
6  LA JAGUA DEL PILAR  44420 MULTIPOLYGON (((-73.08581 1...
7             URUMITA  44855 MULTIPOLYGON (((-73.00927 1...
8          VILLANUEVA  44874 MULTIPOLYGON (((-72.9918 10...
9            RIOHACHA  44001 MULTIPOLYGON (((-72.89025 1...
10            ALBANIA  44035 MULTIPOLYGON (((-72.52248 1...

p.sf <- st_as_sf(precip.points)
(precip.sf = st_intersection(guajira.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 684 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -73.625 ymin: 10.425 xmax: -71.125 ymax: 12.425
geographic CRS: WGS 84
First 10 features:
      MUNIC CODIGO         rain               geometry
13   URIBIA  44847 2.917222e+00 POINT (-71.725 12.425)
13.1 URIBIA  44847 4.168321e-01 POINT (-71.625 12.425)
13.2 URIBIA  44847 6.771785e-15 POINT (-71.575 12.425)
13.3 URIBIA  44847 0.000000e+00 POINT (-71.525 12.425)
13.4 URIBIA  44847 1.536622e+00 POINT (-71.725 12.375)
13.5 URIBIA  44847 6.167189e-02 POINT (-71.675 12.375)
13.6 URIBIA  44847 0.000000e+00 POINT (-71.625 12.375)
13.7 URIBIA  44847 7.484411e-03 POINT (-71.575 12.375)
13.8 URIBIA  44847 2.421079e-02 POINT (-71.525 12.375)
13.9 URIBIA  44847 4.976604e-02 POINT (-71.475 12.375)

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

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

class       : SpatialPointsDataFrame 
features    : 684 
extent      : 1049433, 1321637, 1644788, 1867198  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +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,             rain 
min values  :    ALBANIA,  44001,                0 
max values  : VILLANUEVA,  44874, 177.722457885742 

shapefile(precip2, filename='./chirps/precip2.shp', overwrite=TRUE)
precip2$rainfall <- round(precip2$rain, 1)
precip2

class       : SpatialPointsDataFrame 
features    : 684 
extent      : 1049433, 1321637, 1644788, 1867198  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :      MUNIC, CODIGO,             rain, rainfall 
min values  :    ALBANIA,  44001,                0,        0 
max values  : VILLANUEVA,  44874, 177.722457885742,    177.7 

(guajira2 <- as(guajira.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

class       : SpatialPolygonsDataFrame 
features    : 15 
extent      : 1045090, 1322921, 1641659, 1870870  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +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  :    ALBANIA,  44001 
max values  : VILLANUEVA,  44874 

shapefile(guajira2, filename='./chirps/montes2.shp', overwrite=TRUE)
precip2@bbox <- guajira2@bbox
tm_shape(guajira2) + tm_polygons() +
  tm_shape(precip2) +
  tm_dots(col="rainfall", palette = "RdYlBu", midpoint = 1.0,
             title="Precipitacion muestreada \n(en mm)", size=0.2) +
  tm_text("rainfall", just="center", xmod=.4, size = 0.3) +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in outputCRS object has comment, which is lost in output

library(spatstat.data)
library(nlme)
library(dplyr)
library(raster)
library(rpart)
library(gstat)
library(DirichletReg)
library(Formula)

th <- as(dirichlet(as.ppp(precip2)), "SpatialPolygons")

Error in dirichlet(as.ppp(precip2)) : 
  no se pudo encontrar la función "dirichlet"

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.