Series temporales de imagenes raster. Parte 1
Daniela Ballari
Septiembre de 2018
Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ
Objetivo de la práctica. Explorar la libreria raster a través del uso de imágenes sateliteles multi-temporales de precipitación. Las imágenes a utilizar son de la misión TRMM de la Nasa (http://trmm.gsfc.nasa.gov/3b43.html). Se utilizarán imagenes de precipitación mensual acumulada para el año 2003. El formato de las imágenes es nc. Además se realizará comparación pixel a pixel entre las imágenes y pluviómetros
Temario
A- Cargar librerias y datos
B- Crear raster brick (multicapa)
C- Explorar archivo raster multitemporal (multicapa)
D- Visualizar archivo raster multitemporal (multicapa)
E- Operaciones sobre raster multitemporal (multicapa)
F- Operaciones sobre raster individuales
G- Comparar con estaciones pluviométricas
G.1- Preparar datos de estaciones pluviométricas
G.2- Extraer serie temporal de estaciones e imágenes multitemporales
G.3- Comparar series temporales
A- Cargar librerias y datos
Datos para esta práctica Aqui
rm(list=ls())#borrar todo lo que haya en memoria
#Cargar librarias
library(raster)## Loading required package: sp## Warning: package 'sp' was built under R version 3.4.4library(sf)## Warning: package 'sf' was built under R version 3.4.4## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3library(latticeExtra)## Loading required package: lattice## Loading required package: RColorBrewerlibrary(rasterVis)## Warning: package 'rasterVis' was built under R version 3.4.4library(plyr)# Use "arrange"
library(hydroGOF) # Si hydroGOF no se carga por un error, instalar la libreria "e1071" y volver a cargar hydroGOF## Loading required package: zoo## Warning: package 'zoo' was built under R version 3.4.4##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric#Definir directorio de trabajo, los datos deben estar almacenados en este directorio
setwd("C:/R_ecohidrologia/seriestemporalesraster")
#Limite del área de estudio
ecuador<-st_read("provincias_disuelto.shp") #Fuente SNI, archivo shp## Reading layer `provincias_disuelto' from data source `C:\R_ecohidrologia\seriestemporalesraster\provincias_disuelto.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 490690.4 ymin: 9445216 xmax: 1147852 ymax: 10160820
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defsecuador<-st_transform(ecuador, 4326)
plot(ecuador["DPA_PROVIN"])
B- Crear raster brick (multicapa). Opción 1
#Cargar imagenes anuales y crear un raster brick
files <- list.files(pattern='.G3.nc') #extraer listado de archivos
files## [1] "TRMM_3B43_ACC.007_20030101_000000.G3.nc"
## [2] "TRMM_3B43_ACC.007_20030201_000000.G3.nc"
## [3] "TRMM_3B43_ACC.007_20030301_000000.G3.nc"
## [4] "TRMM_3B43_ACC.007_20030401_000000.G3.nc"
## [5] "TRMM_3B43_ACC.007_20030501_000000.G3.nc"
## [6] "TRMM_3B43_ACC.007_20030601_000000.G3.nc"
## [7] "TRMM_3B43_ACC.007_20030701_000000.G3.nc"
## [8] "TRMM_3B43_ACC.007_20030801_000000.G3.nc"
## [9] "TRMM_3B43_ACC.007_20030901_000000.G3.nc"
## [10] "TRMM_3B43_ACC.007_20031001_000000.G3.nc"
## [11] "TRMM_3B43_ACC.007_20031101_000000.G3.nc"
## [12] "TRMM_3B43_ACC.007_20031201_000000.G3.nc"imageraster <- raster(files[1])#cargar imagen## Loading required namespace: ncdf4spplot(imageraster, col.regions=colorRampPalette(c("white", "blue"))(255))
imageraster <- crop(imageraster, as(ecuador, "Spatial"))#cortar para ecuador
spplot(imageraster, col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador, "Spatial")[1], fill="transparent", col="black", under = FALSE))
imageraster <- mask(imageraster, as(ecuador, "Spatial"))# hacer mascara para ecuador
spplot(imageraster, col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador, "Spatial")[1], fill="transparent", col="black", under = FALSE))
names(imageraster)## [1] "accumulated_r"names(imageraster) <- substr(files[1], 19, 24)#extraer nombre del archivo
names(imageraster)## [1] "X200301"rbrick <- imageraster#asignar a la variable rbrick
#Estructura for
for (i in 1:10)
{
print(i)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10for (i in 2:length(files)) {
imageraster <- raster(files[i])#cargar imagen
imageraster <- crop(imageraster, as(ecuador, "Spatial"))#cortar para ecuador
imageraster <- mask(imageraster, as(ecuador, "Spatial"))# hacer mascara para ecuador
names(imageraster) <- substr(files[i], 19, 24)#extraer nombre del archivo
rbrick <- addLayer(rbrick,imageraster) # acumular las 12 imágenes en una misma variable raster
}
names(rbrick)## [1] "X200301" "X200302" "X200303" "X200304" "X200305" "X200306" "X200307"
## [8] "X200308" "X200309" "X200310" "X200311" "X200312"B- Crear raster brick (multicapa). Opción 2
imageraster <- lapply(files,raster)
spplot(imageraster[[1]], col.regions=colorRampPalette(c("white", "blue"))(255))
imageraster <- lapply(imageraster, crop, y=as(ecuador,"Spatial"))#cortar para ecuador
spplot(imageraster[[1]], col.regions=colorRampPalette(c("white", "blue"))(255)) +
as.layer(spplot(as(ecuador[1], "Spatial"), fill="transparent", col="black", under = FALSE))
C- Explorar archivo raster multitemporal (multicapa)
rbrick## class : RasterStack
## dimensions : 26, 23, 598, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## names : X200301, X200302, X200303, X200304, X200305, X200306, X200307, X200308, X200309, X200310, X200311, X200312
## min values : 0.09069275, 17.36271667, 37.70506668, 6.51883316, 0.07718430, 0.19843662, 0.00000000, 0.00000000, 0.00000000, 0.33404827, 0.00000000, 1.26318693
## max values : 439.1327, 495.3195, 415.3514, 526.6675, 551.2563, 523.4909, 386.9567, 320.7042, 382.8305, 404.6749, 427.7702, 500.5126class(rbrick)## [1] "RasterStack"
## attr(,"package")
## [1] "raster"spplot(rbrick, zlim=c(0,550), col.regions=colorRampPalette(c("white", "blue"))(255))
#Subset
subset(rbrick, 1) #acceder a un raster en particular## class : RasterLayer
## dimensions : 26, 23, 598 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : X200301
## values : 0.09069275, 439.1327 (min, max)rbrick[[1]]## class : RasterLayer
## dimensions : 26, 23, 598 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : X200301
## values : 0.09069275, 439.1327 (min, max)spplot(rbrick[[1]], asp=1, zlim=c(0,550), col.regions=colorRampPalette(c("white", "blue"))(255))
rbrick[20,10]#Serie temporal de precipitación en el pixel 20, 10 (para todas las imágenes)## X200301 X200302 X200303 X200304 X200305 X200306 X200307
## [1,] 88.96648 85.89122 185.123 181.3845 147.9702 174.4233 160.5187
## X200308 X200309 X200310 X200311 X200312
## [1,] 105.3044 144.8949 117.1332 105.5485 100.1129rbrick[20,]# fila 20, todas las columnas## X200301 X200302 X200303 X200304 X200305 X200306
## [1,] NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA
## [4,] 45.57074 62.46648 83.27581 47.65187 2.453969 1.718991
## [5,] 45.77591 82.20290 100.04018 76.68047 21.874407 13.663373
## [6,] 50.57272 90.79436 120.38966 93.53438 27.622101 22.970192
## [7,] 56.63242 93.33072 136.52361 110.78365 40.370876 38.110489
## [8,] 64.43598 91.00589 153.18958 129.77982 63.639210 60.319290
## [9,] 74.81738 84.48792 169.95088 158.80885 110.463890 126.191193
## [10,] 88.96648 85.89122 185.12305 181.38448 147.970184 174.423340
## [11,] 113.87217 120.16452 207.09549 202.61633 193.877625 210.167175
## [12,] NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA
## X200307 X200308 X200309 X200310 X200311 X200312
## [1,] NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA
## [4,] 6.126318 4.303387 7.140658 6.730223 13.72863 50.80203
## [5,] 8.552402 5.711971 19.042469 16.308571 40.99892 51.75758
## [6,] 12.612604 8.162335 26.192497 22.437077 56.37257 60.24499
## [7,] 24.548985 15.806217 45.282528 33.660614 72.16846 72.68435
## [8,] 44.508816 29.165894 77.531677 59.775707 90.41068 85.92162
## [9,] 120.527000 77.938011 125.137634 97.625244 92.57514 84.02791
## [10,] 160.518738 105.304428 144.894928 117.133240 105.54851 100.11290
## [11,] 181.008942 122.201385 178.082489 141.142197 129.76743 127.85788
## [12,] NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA
## [23,] NA NA NA NA NA NArbrick[,10]#Pixel columna 10, todas las filas## X200301 X200302 X200303 X200304 X200305 X200306 X200307
## [1,] NA NA NA NA NA NA NA
## [2,] 173.49321 117.91660 255.37918 515.4109 551.25629 306.53571 144.28146
## [3,] 132.68826 139.93221 181.63205 325.5596 232.40576 141.68549 77.94997
## [4,] 122.68372 152.33156 159.19937 271.7531 143.72093 93.94499 56.98956
## [5,] 117.59333 157.93547 141.01373 234.5410 91.39257 65.43287 45.05561
## [6,] 118.80257 160.70628 132.05896 224.6595 84.54388 63.05284 39.11699
## [7,] 113.25858 139.14919 111.87453 208.6970 71.07278 62.43953 24.60302
## [8,] 124.12768 140.25391 109.28937 203.5284 69.35516 66.25950 20.34694
## [9,] 131.98589 142.90865 110.85727 190.0448 66.58495 71.18283 31.76075
## [10,] 112.16577 131.51364 112.08982 173.7358 62.29975 77.06812 32.16726
## [11,] 99.35641 119.51439 128.50656 157.7056 95.34114 127.10104 74.18602
## [12,] 88.09542 114.10769 132.65137 149.3893 101.96280 145.65305 91.15729
## [13,] 82.92030 106.89030 132.80136 146.5458 109.11610 150.70206 101.89546
## [14,] 74.16940 102.24965 124.66815 141.8571 113.22992 140.35556 87.73564
## [15,] 66.24711 99.69589 101.58868 142.4259 127.35599 117.77525 65.01074
## [16,] 58.19965 95.51315 92.91199 136.6413 132.26224 112.11729 56.00095
## [17,] 55.68230 89.27673 94.28621 138.8637 134.42303 114.93878 61.08647
## [18,] 59.61627 83.52740 110.07909 145.6814 139.20125 124.48808 83.78522
## [19,] 85.91236 93.14344 167.59167 171.7999 144.84970 166.09457 145.61981
## [20,] 88.96648 85.89122 185.12305 181.3845 147.97018 174.42334 160.51874
## [21,] 95.33681 89.91202 199.03651 190.0474 157.90430 178.63654 167.63092
## [22,] 106.83957 107.00260 214.67062 183.0278 163.59276 187.75493 168.21893
## [23,] 125.54922 139.10556 218.47600 188.8070 196.83110 199.40977 155.60075
## [24,] NA NA NA NA NA NA NA
## [25,] NA NA NA NA NA NA NA
## [26,] NA NA NA NA NA NA NA
## X200308 X200309 X200310 X200311 X200312
## [1,] NA NA NA NA NA
## [2,] 258.32419 277.79495 283.35669 406.13901 342.69077
## [3,] 102.18410 150.83650 197.36362 281.60498 205.23729
## [4,] 60.31007 107.37806 149.24821 218.77335 154.64313
## [5,] 37.06458 79.37399 118.44450 175.99127 127.60810
## [6,] 32.58199 73.31939 108.40964 160.63510 120.01450
## [7,] 23.63161 67.84982 98.94843 135.95518 98.55187
## [8,] 24.02166 65.39985 94.18610 129.98277 93.31880
## [9,] 28.62406 70.18970 93.58798 117.33698 98.29994
## [10,] 27.25663 71.62355 87.43563 113.42832 102.62121
## [11,] 53.98546 82.89888 94.52548 102.33572 117.69818
## [12,] 62.41408 88.04301 98.44228 99.39512 122.16456
## [13,] 68.38576 97.28098 96.87968 95.28327 127.16487
## [14,] 62.25468 93.96672 89.59691 89.26069 120.90095
## [15,] 52.05892 82.84716 72.00885 79.34607 120.85791
## [16,] 47.87391 83.63474 66.42889 71.07800 116.15646
## [17,] 52.32886 88.16235 66.33945 74.75337 111.88361
## [18,] 59.36547 99.97610 77.65502 82.88043 110.97887
## [19,] 91.66550 139.97279 106.39729 99.27980 104.77837
## [20,] 105.30443 144.89493 117.13324 105.54851 100.11290
## [21,] 110.53720 156.85950 122.89774 112.03171 99.21403
## [22,] 110.37660 152.02502 128.55667 121.32169 103.75388
## [23,] 107.83009 156.29993 133.91838 141.45461 124.59822
## [24,] NA NA NA NA NA
## [25,] NA NA NA NA NA
## [26,] NA NA NA NA NArbrick[[1]][20,10]#un determinado pixel de un raster en particular##
## 88.96648D- Visualizacion de archivo raster multitemporal (multicapa)
levelplot(rbrick[[1]]) #mean # funcion de rasterVis
#color adecuado a precipitación
mytheme <- rasterTheme(region = brewer.pal(9, "Blues"))
levelplot(rbrick[[1]], par.settings = mytheme)
levelplot(rbrick[[1]], par.settings = mytheme, margin=F)#eliminar los margenes
levelplot(rbrick[[1]], par.settings = mytheme, margin = list(FUN = 'sd'))#plotar en margenes otra funcion
levelplot(rbrick[[1]], par.settings = mytheme,
margin = list(FUN = 'mean'),contour = TRUE) # Isolineas y promedio en margenes
levelplot(rbrick[[1]], par.settings = mytheme, margin = list(FUN = 'mean'),
contour = TRUE)+ layer(sp.polygons(as(ecuador,"Spatial"), lwd=2, col='black')) # Plotar con poligono de ecuador, + layer es funcion de libreria lattice
#plotear imagenes multitemporales
levelplot(rbrick[[1:4]])
levelplot(rbrick[[1:3]], par.settings = mytheme,colorkey = list(space = "bottom"),
layout = c(3, 1, 1)) + layer(as(ecuador,"Spatial"))
E- Operaciones sobre raster multitemporal (multicapa)
#Precipitación anual acumulada
precipitacion <- calc(rbrick, fun=sum,na.rm=TRUE)
precipitacion## class : RasterLayer
## dimensions : 26, 23, 598 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0, 4856.859 (min, max)names(precipitacion) <- "sum"
#Precipitación media anual
precipitacion$media <- calc(rbrick, fun=mean,na.rm=TRUE)
#Desviación estándar
precipitacion$ds <- calc(rbrick, fun=sd,na.rm=TRUE)
precipitacion## class : RasterStack
## dimensions : 26, 23, 598, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## names : sum, media, ds
## min values : 0.00000, 12.48985, 25.04007
## max values : 4856.8588, 404.7382, 186.5787spplot(precipitacion, asp=1, col.regions=colorRampPalette(c("white", "blue"))(255))
spplot(precipitacion[[1]], asp=1, col.regions=colorRampPalette(c("white", "blue"))(255), main="Suma")
spplot(precipitacion[[2]], asp=1, col.regions=colorRampPalette(c("white", "blue"))(255), main="Media")
spplot(precipitacion[[3]], asp=1, col.regions=colorRampPalette(c("white", "blue"))(255), main="Desviación estándar")
#Correlación (pearson) en raster multicapas
layerStats(rbrick, 'pearson',na.rm=TRUE)## $`pearson correlation coefficient`
## X200301 X200302 X200303 X200304 X200305 X200306
## X200301 1.0000000 0.83848180 0.6154167 0.7911187 0.4439985 0.2310618
## X200302 0.8384818 1.00000000 0.4993102 0.6370454 0.3989530 0.1830643
## X200303 0.6154167 0.49931017 1.0000000 0.8723458 0.8087434 0.7385603
## X200304 0.7911187 0.63704535 0.8723458 1.0000000 0.8345786 0.6844630
## X200305 0.4439985 0.39895299 0.8087434 0.8345786 1.0000000 0.9335059
## X200306 0.2310618 0.18306431 0.7385603 0.6844630 0.9335059 1.0000000
## X200307 0.2428997 0.14264062 0.7330597 0.6756128 0.8894881 0.9592494
## X200308 0.2181924 0.16039647 0.7523022 0.7088524 0.9326963 0.9553273
## X200309 0.2249863 0.04453638 0.6673681 0.6751871 0.8544857 0.8926644
## X200310 0.2421112 0.18418498 0.7545754 0.7348869 0.9212326 0.9424297
## X200311 0.2028289 0.03485573 0.6843069 0.6878389 0.8403926 0.8979400
## X200312 0.4559963 0.27925977 0.7786102 0.8159339 0.9165782 0.8994396
## X200307 X200308 X200309 X200310 X200311 X200312
## X200301 0.2428997 0.2181924 0.22498628 0.2421112 0.20282892 0.4559963
## X200302 0.1426406 0.1603965 0.04453638 0.1841850 0.03485573 0.2792598
## X200303 0.7330597 0.7523022 0.66736814 0.7545754 0.68430688 0.7786102
## X200304 0.6756128 0.7088524 0.67518710 0.7348869 0.68783892 0.8159339
## X200305 0.8894881 0.9326963 0.85448566 0.9212326 0.84039263 0.9165782
## X200306 0.9592494 0.9553273 0.89266440 0.9424297 0.89793996 0.8994396
## X200307 1.0000000 0.9337476 0.93849358 0.9045894 0.87914645 0.8984797
## X200308 0.9337476 1.0000000 0.91107327 0.9612322 0.91296105 0.8733623
## X200309 0.9384936 0.9110733 1.00000000 0.8758349 0.90912844 0.9262050
## X200310 0.9045894 0.9612322 0.87583489 1.0000000 0.92675990 0.8640036
## X200311 0.8791465 0.9129610 0.90912844 0.9267599 1.00000000 0.8781094
## X200312 0.8984797 0.8733623 0.92620498 0.8640036 0.87810945 1.0000000
##
## $mean
## X200301 X200302 X200303 X200304 X200305 X200306 X200307
## 157.70113 210.64887 201.39013 248.86259 224.71408 185.75193 107.97345
## X200308 X200309 X200310 X200311 X200312
## 92.39613 114.35688 140.42800 150.62266 177.51512pairs(rbrick)#ploteo de correlaciones
F- Operaciones sobre raster individuales
nlayers(rbrick)## [1] 12summary(rbrick[[1]])## X200301
## Min. 0.09069275
## 1st Qu. 86.96751022
## Median 132.16464233
## 3rd Qu. 209.91656494
## Max. 439.13269043
## NA's 275.00000000rbrick[[1]]## class : RasterLayer
## dimensions : 26, 23, 598 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : X200301
## values : 0.09069275, 439.1327 (min, max)spplot(rbrick[[1]], asp=1, col.regions=colorRampPalette(c("white", "blue"))(255), main="Enero")
Moran(rbrick[[1]])#autocorrelación espacial## [1] 0.8019856Moran(rbrick[[3]])## [1] 0.7580466moran.ts<- sapply(imageraster, Moran)
barplot(moran.ts)
hist(rbrick[[1]])#histograma
boxplot(rbrick[[1]])#grafico de caja
plot(rbrick[[1]])
# click(rbrick[[1]])#click en Esc para salir de la opción
#Estadística descriptiva
cellStats(rbrick[[1]], stat="sum",na.rm=TRUE)## [1] 50937.46#Funcion para agregar una fila mensual con media, acumulado y desviación estandar
#Estructura del data.frame
precipitacion_mensual <- data.frame(suma=numeric(), media=numeric(), ds=numeric())
for (i in 1:nlayers(rbrick))
{
precipitacion_mensual<- rbind(precipitacion_mensual,
c(cellStats(rbrick[[i]], stat="sum",na.rm=TRUE),
cellStats(rbrick[[i]], stat="mean",na.rm=TRUE) ,
cellStats(rbrick[[i]], stat="sd",na.rm=TRUE) ))
}
names(precipitacion_mensual) <- c("suma", "media", "ds")
precipitacion_mensual## suma media ds
## 1 50937.46 157.70113 94.11612
## 2 68039.58 210.64887 94.89979
## 3 65049.01 201.39013 83.03529
## 4 80382.62 248.86259 122.17373
## 5 72582.65 224.71408 164.55182
## 6 59997.87 185.75193 167.37007
## 7 34875.43 107.97345 97.81509
## 8 29843.95 92.39613 87.52096
## 9 36937.27 114.35688 98.09942
## 10 45358.24 140.42800 109.89811
## 11 48651.12 150.62266 120.33863
## 12 57337.38 177.51512 122.07235prec.mensual.sum<- sapply(imageraster, cellStats,stat="sum",na.rm=TRUE)
prec.mensual.mean<- sapply(imageraster, cellStats,stat="mean",na.rm=TRUE)
prec.mensual.sd <- sapply(imageraster, cellStats,stat="sd",na.rm=TRUE)
cbind(prec.mensual.sum,prec.mensual.mean,prec.mensual.sd)## prec.mensual.sum prec.mensual.mean prec.mensual.sd
## [1,] 93016.48 155.5460 90.41711
## [2,] 122342.69 204.5864 91.83003
## [3,] 120871.79 202.1267 84.87592
## [4,] 140575.05 235.0753 126.85509
## [5,] 143592.00 240.1204 165.79499
## [6,] 117536.56 196.5494 155.45141
## [7,] 69469.60 116.1699 93.32025
## [8,] 63421.96 106.0568 89.57575
## [9,] 78073.82 130.5582 102.99934
## [10,] 88307.68 147.6717 111.08084
## [11,] 93405.02 156.1957 115.62527
## [12,] 115353.78 192.8993 128.42079#Indice de autocorrelacion local LISA.
ml<- MoranLocal(rbrick[[1]])
ml## class : RasterLayer
## dimensions : 26, 23, 598 (nrow, ncol, ncell)
## resolution : 0.25, 0.25 (x, y)
## extent : -81, -75.25, -5, 1.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : -0.08601961, 7.619742 (min, max)levelplot(ml, margin=FALSE)
#correlación local, pixel a pixel
for (i in 1:ncell(rbrick))
{
cor<- data.frame(cor=Moran(raster(rbrick[i])), xyFromCell(rbrick, i))
if(i==1) (cor.accum = cor)
if(i!=1) (cor.accum = rbind(cor.accum, cor))
}
names(cor.accum)## [1] "cor" "x" "y"coordinates(cor.accum)<- c("x","y")
gridded(cor.accum)<- TRUE
plot(raster(cor.accum))
# Correlacion entre pixels
rbrick.df<- as(rbrick, "SpatialPointsDataFrame")
rbrick.df## class : SpatialPointsDataFrame
## features : 323
## extent : -80.875, -75.375, -4.875, 1.125 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## variables : 12
## names : X200301, X200302, X200303, X200304, X200305, X200306, X200307, X200308, X200309, X200310, X200311, X200312
## min values : 0.0906927511096001, 17.3627166748047, 37.7050666809082, 6.51883316040039, 0.0771843045949936, 0.198436617851257, 0, 0, 0, 0.334048271179199, 0, 1.26318693161011
## max values : 439.132690429688, 495.319458007812, 415.351409912109, 526.66748046875, 551.256286621094, 523.490905761719, 386.956665039062, 320.704223632812, 382.830474853516, 404.674896240234, 427.770172119141, 500.512573242188head(rbrick.df)## X200301 X200302 X200303 X200304 X200305 X200306 X200307
## 1 159.4487 129.9251 209.4438 319.6915 266.3918 129.74361 65.73835
## 2 180.7647 128.9594 251.6386 466.7869 470.2444 256.05817 120.19913
## 3 173.4932 117.9166 255.3792 515.4109 551.2563 306.53571 144.28146
## 4 157.6488 113.6176 241.5055 496.6528 546.6296 314.25665 145.34225
## 5 301.7777 272.6137 256.0273 298.6614 224.0011 72.95753 42.74660
## 6 290.4550 255.1562 267.1194 328.2341 264.7957 92.93478 49.59270
## X200308 X200309 X200310 X200311 X200312
## 1 125.87242 101.53854 108.17146 150.92163 155.3847
## 2 225.47449 224.79810 219.78867 323.55023 293.6209
## 3 258.32419 277.79495 283.35669 406.13901 342.6908
## 4 246.75623 287.26251 302.73969 427.77017 350.1609
## 5 59.97577 46.15705 77.56047 48.33755 142.8242
## 6 72.28557 61.71827 97.86996 78.91939 162.5535cor(t(rbrick.df@data[50,]), t(rbrick.df@data[60,]))## 60
## 50 0.6291161cor(t(rbrick.df@data[50,]), t(rbrick.df@data[60,]), method="spearman")## 60
## 50 0.5104895G- Comparar con estaciones pluviométricas
G.1- Preparar datos de estaciones pluviométricas
estaciones<-st_read("estaciones_32717.shp") #estaciones metereológicas## Reading layer `estaciones_32717' from data source `C:\R_ecohidrologia\seriestemporalesraster\estaciones_32717.shp' using driver `ESRI Shapefile'
## Simple feature collection with 40 features and 14 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 586178.3 ymin: 9561300 xmax: 731094.3 ymax: 9667788
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defsestaciones2<- estaciones[,c("CODIGO","NOMBRE_DE","Lat","Long", "x","y")]
names(estaciones2)<- c("estacion","nombre", "lat", "long", "x", "y") #asignar nuevos nombres a las columnas
zona_estudio<-st_read("prov_zona_estudio.shp")## Reading layer `prov_zona_estudio' from data source `C:\R_ecohidrologia\seriestemporalesraster\prov_zona_estudio.shp' using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 557012.9 ymin: 9445216 xmax: 979220 ymax: 9841834
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs### Cargue un archivo de texto con observaciones y elimine valores vacíos
observaciones <- read.table("observaciones_2003.csv")
table(complete.cases(observaciones$VALOR))##
## FALSE TRUE
## 1786 6244observaciones2 <- observaciones[complete.cases(observaciones$VALOR),]#excluir valores vacios
table(complete.cases(observaciones2$VALOR))##
## TRUE
## 6244Asigne coordenadas a las observaciones. Proceda a unir los data.frame estaciones2 y observaciones2 con la función mergepara que las observaciones hereden las coordenadas de sus respectivas estaciones. Las columnas a unir son observaciones2$ESTACION y estaciones2$estacion.
observ_ref <- merge(observaciones2, estaciones2,by.x = "ESTACION", by.y = "estacion")La variable obtenida es de tipo data.frame, es decir no está georreferenciada. Sin embargo, ahora el data.frame contiene una columna llamada x e y que se utilizarán para georreferenciar con la ayuda de la función coordinates. Luego, se asignará un CRS ya que el SpatialPointsDataFrame no tiene asignado ninguno.
Aregar observaciones mensuales
head(observ_ref)## ESTACION FECHA VALOR STATUS year day month nombre lat
## 1 M012 01/01/2003 0 2003 1 1 LA CUCA -3.492222
## 2 M012 02/01/2003 0 2003 2 1 LA CUCA -3.492222
## 3 M012 03/01/2003 0 2003 3 1 LA CUCA -3.492222
## 4 M012 04/01/2003 0 2003 4 1 LA CUCA -3.492222
## 5 M012 05/01/2003 0 2003 5 1 LA CUCA -3.492222
## 6 M012 06/01/2003 0 2003 6 1 LA CUCA -3.492222
## long x y NA
## 1 -80.05694 604749.6 9613947 POINT (604749.6 9613947)
## 2 -80.05694 604749.6 9613947 POINT (604749.6 9613947)
## 3 -80.05694 604749.6 9613947 POINT (604749.6 9613947)
## 4 -80.05694 604749.6 9613947 POINT (604749.6 9613947)
## 5 -80.05694 604749.6 9613947 POINT (604749.6 9613947)
## 6 -80.05694 604749.6 9613947 POINT (604749.6 9613947)p.2003.mensual <- aggregate(observ_ref$VALOR,by=list(observ_ref$ESTACION,observ_ref$month), FUN=sum)
p.2003.unique<- unique(observ_ref[,c("ESTACION", "x", "y")])#estaciones unicas
head(p.2003.mensual)## Group.1 Group.2 x
## 1 M012 1 42.7
## 2 M040 1 32.0
## 3 M142 1 26.4
## 4 M180 1 97.3
## 5 M185 1 6.8
## 6 M232 1 151.0head(p.2003.unique)## ESTACION x y
## 1 M012 604749.6 9613947
## 353 M040 635322.0 9631878
## 718 M142 696350.5 9599623
## 1069 M180 653688.8 9591251
## 1413 M185 640761.4 9662797
## 1539 M232 602296.0 9568778p.2003.mensual.unido<- merge(x = p.2003.mensual, y = p.2003.unique, by.x = "Group.1", by.y="ESTACION", all.x=TRUE)
coordinates(p.2003.mensual.unido) <- c("x.y","y")
utm17 <- "+init=epsg:32717"
proj4string(p.2003.mensual.unido) <- CRS(utm17)
names(p.2003.mensual.unido)<- c("ESTACION", "MES", "PREC")
head(p.2003.mensual.unido)## ESTACION MES PREC
## 1 M012 1 42.7
## 2 M012 10 0.0
## 3 M012 4 0.0
## 4 M012 8 0.0
## 5 M012 11 0.0
## 6 M012 7 9.4bubble(p.2003.mensual.unido, "PREC", col="black", do.sqrt=FALSE, pch=21,
main="PREC mensual 2003", sp.layout=list("sp.polygons",
col="grey", fill="transparent", as(zona_estudio, "Spatial")))
G.2- Extraer serie temporal de estaciones e imágenes multitemporales
#Extraer serie temporal de estaciones
count=0
for (i in unique(p.2003.mensual.unido$ESTACION)) #
{
count=count + 1
data <-data.frame(subset(p.2003.mensual.unido,p.2003.mensual.unido$ESTACION == i))
row<- data.frame(i, t(arrange(data, MES)[,3] ) )#arrange (plyr package)-> orderar datos por columnas
coord <- data.frame(i, unique(data)[1,c(4,5)]) #extraer coordenada de la serie temporal
if(length(row)!=13) next #para omitir estaciones con meses sin datos
if(count ==1) (time.serie = row)
if(count ==1) (coordacum <- coord)
if(count !=1) (time.serie = rbind(time.serie, row))
if(count !=1) (coordacum <- rbind(coordacum, coord))
}
time.serie.est<- time.serie
time.serie.est## i X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11
## 1 M012 42.7 72.1 18.0 0.0 6.2 36.0 9.4 0.0 0.0 0.0 0.0
## 2 M040 32.0 73.7 110.4 10.3 39.9 19.8 39.1 30.3 24.8 48.8 40.8
## 3 M142 26.4 48.6 119.5 97.4 54.0 28.9 4.7 13.6 46.8 43.8 65.3
## 4 M180 97.3 90.5 264.9 199.3 53.3 15.1 4.7 8.2 0.6 17.4 92.8
## 5 M232 151.0 198.6 246.8 95.9 23.7 5.1 1.0 0.0 0.0 0.0 14.1
## 6 M239 80.1 128.2 242.3 133.5 21.2 3.8 3.4 7.2 0.0 10.0 60.4
## 7 M292 34.8 69.8 76.3 3.7 7.1 18.2 8.8 14.2 20.8 20.6 27.4
## 8 M418 18.9 51.8 93.4 64.1 35.3 75.4 30.3 6.8 37.8 28.2 83.4
## 9 M420 10.4 32.8 118.3 83.9 26.4 51.7 16.2 2.2 32.2 46.6 64.9
## 10 M421 7.2 0.0 60.1 23.4 33.8 10.5 20.4 11.3 23.2 1.3 39.5
## 11 M422 28.4 29.5 65.2 62.3 9.9 18.4 13.5 0.0 1.4 11.6 13.4
## 12 M481 73.2 71.0 96.4 29.9 32.0 20.4 10.2 15.2 7.3 21.0 37.5
## 13 M503 134.7 110.8 143.4 252.6 517.4 252.8 217.9 124.4 142.6 118.9 85.3
## 14 M760 127.3 174.4 277.5 264.3 26.8 18.5 3.0 0.0 0.0 7.2 6.3
## 15 M773 148.6 168.5 220.8 154.3 39.4 20.7 15.5 7.2 1.5 5.0 44.7
## X12
## 1 39.8
## 2 22.1
## 3 79.4
## 4 96.0
## 5 104.1
## 6 54.7
## 7 28.2
## 8 100.1
## 9 39.9
## 10 5.6
## 11 15.2
## 12 65.5
## 13 199.2
## 14 96.8
## 15 97.4#Extraer serie temporal de raster brick
names(coordacum)## [1] "i" "x.y" "y"coordinates(coordacum) <- c("x.y","y")
coordacum_spdf <- as(coordacum, "SpatialPoints")
projection(rbrick) <- CRS("+init=epsg:4326")
rbrick_utm17 <- projectRaster(rbrick,crs=utm17)#transformación utilizando libreria raster
#Extraer serie temporal en base a coordenadas
time.serie.raster <- extract(rbrick_utm17, coordacum_spdf)
time.serie.raster <- data.frame(coordacum@data[,1], time.serie.raster)
time.serie.raster## coordacum.data...1. X200301 X200302 X200303 X200304 X200305
## 1 M012 45.44496 62.43122 83.48297 47.72526 2.522950
## 2 M040 45.70163 82.09046 100.32286 76.75239 21.922318
## 3 M142 62.99831 83.00013 167.08805 131.78493 66.716921
## 4 M180 44.47901 79.61869 141.44620 95.55386 27.315649
## 5 M232 28.72612 61.19318 110.14855 46.32382 1.721534
## 6 M239 44.47901 79.61869 141.44620 95.55386 27.315649
## 7 M292 45.70163 82.09046 100.32286 76.75239 21.922318
## 8 M418 66.93277 111.30080 143.44410 138.66898 65.838135
## 9 M420 64.48962 90.88616 153.43065 129.98691 63.971696
## 10 M421 64.48962 90.88616 153.43065 129.98691 63.971696
## 11 M422 56.62696 93.17307 136.82355 110.91223 40.507525
## 12 M481 50.52602 90.66223 120.74497 93.65086 27.684701
## 13 M503 60.26559 85.20531 179.89486 134.02856 70.391783
## 14 M760 34.38742 71.10897 135.71106 72.56515 19.352215
## 15 M773 44.47901 79.61869 141.44620 95.55386 27.315649
## X200306 X200307 X200308 X200309 X200310 X200311 X200312
## 1 1.783009 6.137465 4.306101 7.181476 6.770403 13.83710 50.84075
## 2 13.755419 8.588589 5.733906 19.074717 16.352314 41.11461 51.83860
## 3 65.024740 47.842044 30.486745 88.036947 66.718569 95.49937 86.08550
## 4 24.222106 16.184303 10.144355 28.740364 27.072823 64.73616 65.88040
## 5 3.151167 6.389838 3.529628 6.218371 6.621002 16.40956 49.47577
## 6 24.222106 16.184303 10.144355 28.740364 27.072823 64.73616 65.88040
## 7 13.755419 8.588589 5.733906 19.074717 16.352314 41.11461 51.83860
## 8 55.529607 38.857287 25.875161 63.865893 49.862408 79.67368 85.72052
## 9 60.790478 45.033199 29.494535 77.942596 60.087890 90.47382 85.91062
## 10 60.790478 45.033199 29.494535 77.942596 60.087890 90.47382 85.91062
## 11 38.237448 24.680319 15.904433 45.541448 33.893443 72.38210 72.82530
## 12 23.065239 12.721603 8.228354 26.325408 22.555908 56.56415 60.38289
## 13 71.348619 47.996989 29.847855 77.231174 66.278508 94.73285 88.58679
## 14 18.306837 10.492302 4.875886 16.599913 16.393140 42.58568 53.69781
## 15 24.222106 16.184303 10.144355 28.740364 27.072823 64.73616 65.88040F.3- Comparar series temporales
#Calcular bias, RMSE y coeficiente de Pearson's
#Analisis temporal
#library(hydroGOF)
RMSE<-rmse(time.serie.raster[,2:13], time.serie.est[,2:13])
#rmse = sqrt( mean( (sim - obs)^2, na.rm = TRUE) )
plot(RMSE,type="p" )
lines(RMSE)
bias <- pbias(time.serie.raster[,2:13],time.serie.est[,2:13])
#PBIAS = 100 * [ sum( sim - obs ) / sum( obs ) ]
# Percent bias (PBIAS) measures the average tendency of the simulated values to be larger or smaller than their observed ones.
# The optimal value of PBIAS is 0.0, with low-magnitude values indicating accurate model simulation. Positive values indicate overestimation bias, whereas negative values indicate model underestimation bias
plot(bias, type="p")
lines(bias)
gof <- gof(time.serie.raster[,2:13],time.serie.est[,2:13] )## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !## Warning: 'rNSE' can not be computed: some elements in 'obs' are zero !## Warning: 'rd' can not be computed: some elements in 'obs' are zero !gof## X200301 X200302 X200303 X200304 X200305 X200306 X200307 X200308
## ME -16.88 -5.16 -9.61 0.06 -25.20 -6.47 -3.15 -1.11
## MAE 50.43 51.41 69.35 74.77 46.87 26.62 24.04 16.83
## MSE 3700.24 4113.13 6481.47 7174.58 13742.77 2625.36 2246.24 781.19
## RMSE 60.83 64.13 80.51 84.70 117.23 51.24 47.39 27.95
## NRMSE % 116.50 111.30 95.50 97.70 92.40 83.00 87.80 90.00
## PBIAS % -25.00 -5.90 -6.70 0.10 -40.80 -16.30 -11.90 -6.90
## RSR 1.17 1.11 0.96 0.98 0.92 0.83 0.88 0.90
## rSD 0.22 0.22 0.31 0.35 0.19 0.38 0.31 0.35
## NSE -0.46 -0.33 0.02 -0.02 0.09 0.26 0.17 0.13
## mNSE -0.11 -0.12 0.03 -0.09 0.23 0.23 0.13 -0.03
## rNSE -11.27 -Inf -2.88 -Inf 0.66 -0.67 -1.53 -Inf
## d 0.18 0.05 0.36 0.36 0.27 0.55 0.44 0.43
## md 0.19 0.08 0.24 0.20 0.47 0.53 0.44 0.36
## rd -5.92 -Inf -1.55 -Inf 0.73 -0.02 -0.72 -Inf
## cp -1.50 -0.90 -0.08 0.12 0.57 0.67 0.64 0.60
## r -0.67 -0.63 0.22 0.14 0.44 0.55 0.44 0.37
## R2 0.44 0.40 0.05 0.02 0.19 0.30 0.20 0.13
## bR2 0.19 0.25 0.03 0.01 0.04 0.12 0.06 0.04
## KGE -0.86 -0.81 -0.04 -0.07 -0.07 0.22 0.10 0.09
## VE 0.25 0.42 0.52 0.24 0.24 0.33 0.09 -0.05
## X200309 X200310 X200311 X200312
## ME 18.15 8.19 16.88 -1.55
## MAE 27.86 20.10 21.12 37.05
## MSE 1120.94 667.88 734.60 2188.85
## RMSE 33.48 25.84 27.10 46.79
## NRMSE % 90.90 84.00 91.40 94.60
## PBIAS % 80.30 32.30 37.50 -2.20
## RSR 0.91 0.84 0.91 0.95
## rSD 0.77 0.70 0.90 0.30
## NSE 0.11 0.24 0.11 0.04
## mNSE -0.22 0.06 0.13 0.03
## rNSE -Inf -Inf -Inf -30.82
## d 0.71 0.72 0.76 0.39
## md 0.41 0.48 0.58 0.28
## rd -Inf -Inf -Inf -19.11
## cp 0.61 0.65 0.61 0.32
## r 0.63 0.58 0.70 0.22
## R2 0.40 0.33 0.49 0.05
## bR2 0.34 0.26 0.42 0.03
## KGE 0.09 0.39 0.51 -0.04
## VE -0.23 0.21 0.53 0.47gof["r",]## X200301 X200302 X200303 X200304 X200305 X200306 X200307 X200308 X200309
## -0.67 -0.63 0.22 0.14 0.44 0.55 0.44 0.37 0.63
## X200310 X200311 X200312
## 0.58 0.70 0.22plot(gof["r",], type="p")
lines(gof["r",])
# ggof(time.serie.raster[,2:13],time.serie.est[,2:13] )
#Analisis por estaciones
#Estación
i=1
time.serie.est[i,1]## [1] M012
## 15 Levels: M012 M040 M142 M180 M232 M239 M292 M418 M420 M421 M422 ... M773RMSE<-rmse(time.serie.raster[i,2:13], time.serie.est[i,2:13], na.rm=TRUE)
plot(RMSE,type="p", main=paste("estacion", time.serie.est[i,1]) )
lines(RMSE)
bias <- pbias(time.serie.raster[i,2:13],time.serie.est[i,2:13], na.rm=TRUE)## Warning in pbias.default(x[, i], y[, i], na.rm = na.rm, ...):
## 'sum((obs)=0', it is not possible to compute 'pbias'## Warning in pbias.default(x[, i], y[, i], na.rm = na.rm, ...):
## 'sum((obs)=0', it is not possible to compute 'pbias'
## Warning in pbias.default(x[, i], y[, i], na.rm = na.rm, ...):
## 'sum((obs)=0', it is not possible to compute 'pbias'
## Warning in pbias.default(x[, i], y[, i], na.rm = na.rm, ...):
## 'sum((obs)=0', it is not possible to compute 'pbias'
## Warning in pbias.default(x[, i], y[, i], na.rm = na.rm, ...):
## 'sum((obs)=0', it is not possible to compute 'pbias'plot(bias, type="p", main=paste("estacion", time.serie.est[i,1]))
lines(bias)
# gof <- gof(time.serie.raster[i,2:13],time.serie.est[i,2:13], na.rm=TRUE )
# gof
# gof["r",]
# plot(gof["r",], type="p",main=paste("estacion", time.serie.est[i,1]))
# lines(gof["r",])