Analisis_espacial.R
usuario
Thu May 25 08:50:39 2017
#==========================================================================@title Introduccion a analisis de datos espacio-temporales en R @author Ing. Carlos A. Fernandez Palomino @Email cafpxl@gmail.com @details Uso Libre @description Este script comprende varios topicos:
#==========================================================================
# ------------------------------------------------------------------------
# Instalar y cargar paquetes ----------------------------------------
# ------------------------------------------------------------------------
#install.packages("sp")
library(easypackages)
libraries("sp", "rgdal", "raster","lattice","latticeExtra")## Loading required package: sp## Loading required package: rgdal## rgdal: version: 1.2-7, (SVN revision 660)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Program Files/R/R-3.4.0/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Program Files/R/R-3.4.0/library/rgdal/proj
## Linking to sp version: 1.2-4## Loading required package: raster## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewer## All packages loaded successfully# configurar directorio de trabajo
setwd("E:/SENAMHI_2017/Curso_R_SENAMHI/3_Datos_espaciales")
# ------------------------------------------------------------------------
# Analisis de datos espaciales ----------------------------------------
# ------------------------------------------------------------------------
library(sp)
getClass('Spatial')## Class "Spatial" [package "sp"]
##
## Slots:
##
## Name: bbox proj4string
## Class: matrix CRS
##
## Known Subclasses:
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2#==== Crear objetos SpatialPoints y SpatialPointsDataFrame
# Desde cero
xy <- matrix(data=runif(n=100), ncol=2)
head(xy)## [,1] [,2]
## [1,] 0.8246055 0.48237894
## [2,] 0.1405399 0.11955092
## [3,] 0.9977758 0.58527534
## [4,] 0.8677272 0.09923827
## [5,] 0.4737037 0.11361835
## [6,] 0.7188797 0.21040632xySp <- SpatialPoints(xy)
summary(xySp)## Object of class SpatialPoints
## Coordinates:
## min max
## coords.x1 0.007501609 0.9977758
## coords.x2 0.025719077 0.9975398
## Is projected: NA
## proj4string : [NA]
## Number of points: 50plot(xySp, axes=T, asp=1, pch=16)
# Desde un dataframe con columnas x e y
# opcion 1
data <- read.csv("Dato_espacial.csv")
head(data,2)## x y Enero Febrero Marzo Abril Mayo Junio Julio Agosto
## 1 279039.0 8663781 0.33 0.59 0.15 0.14 0.65 1.03 1.25 2.46
## 2 256723.4 8731110 1.07 2.77 2.15 0.39 0.92 3.61 3.34 5.66
## Septiembre Octubre Noviembre Diciembre DEF MAM JJA SON Anual
## 1 2.12 0.78 0.64 0.93 1.85 0.94 4.74 3.54 11.07
## 2 2.23 0.44 0.43 1.40 5.24 3.46 12.61 3.10 24.41
## est id Z
## 1 MODELO.CAMPO.DE.MART 1 131
## 2 DONOSO 2 158xySpD <- SpatialPointsDataFrame(data[,c(1,2)], data)
plot(xySpD, axes=T, asp=1, pch=16)
# opcion 2
data <- read.csv("Dato_espacial.csv")
coordinates(data) <- ~x + y #alt 126 para "~"
class(data)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"plot(data, axes=T, asp=1, pch=16)
# para acceder a la tabla de atributos, utilizar 'data slot' (con el caracter '@' )
nrow(data@data)## [1] 36head(data@data) # que paso con las coordenadas?## Enero Febrero Marzo Abril Mayo Junio Julio Agosto Septiembre Octubre
## 1 0.33 0.59 0.15 0.14 0.65 1.03 1.25 2.46 2.12 0.78
## 2 1.07 2.77 2.15 0.39 0.92 3.61 3.34 5.66 2.23 0.44
## 3 0.66 1.09 0.59 0.67 1.05 1.52 1.82 2.14 1.50 0.69
## 4 1.58 3.27 1.07 0.22 0.56 1.31 2.12 1.88 0.88 0.27
## 5 5.97 11.69 5.89 0.45 0.16 0.01 0.00 0.07 0.06 0.18
## 6 4.41 8.04 4.57 0.77 0.21 0.00 0.02 0.01 0.07 0.05
## Noviembre Diciembre DEF MAM JJA SON Anual est id
## 1 0.64 0.93 1.85 0.94 4.74 3.54 11.07 MODELO.CAMPO.DE.MART 1
## 2 0.43 1.40 5.24 3.46 12.61 3.10 24.41 DONOSO 2
## 3 0.77 0.65 2.40 2.31 5.48 2.96 13.15 VON.HUMBOLDT 3
## 4 0.56 1.39 6.24 1.85 5.31 1.71 15.11 HUAYAN 4
## 5 0.37 2.78 20.44 6.50 0.08 0.61 27.63 SANTA.EULALIA 5
## 6 0.37 1.06 13.51 5.55 0.03 0.49 19.58 CHOSICA 6
## Z
## 1 131
## 2 158
## 3 234
## 4 351
## 5 959
## 6 1220head(coordinates(data),2)# utilizar la funcion 'coordinates' para recuperar las coordenadas## x y
## 1 279039.0 8663781
## 2 256723.4 8731110# plot de estaciones
plot(data, axes=T, asp=1, pch=16)
# filtros sobre los atributos
# Seleccionar solo las estaciones sobre los 500 msnm
est.cuenca.alta <- data[data@data$Z > 500, ]
plot(data, axes=T, asp=1, pch=16)
plot(est.cuenca.alta, pch=1, col="red", cex=3, add=TRUE)
# agregar un nuevo atributo a la columna de atributos
data@data[["Pais"]] <- c("peru")
head(data@data)## Enero Febrero Marzo Abril Mayo Junio Julio Agosto Septiembre Octubre
## 1 0.33 0.59 0.15 0.14 0.65 1.03 1.25 2.46 2.12 0.78
## 2 1.07 2.77 2.15 0.39 0.92 3.61 3.34 5.66 2.23 0.44
## 3 0.66 1.09 0.59 0.67 1.05 1.52 1.82 2.14 1.50 0.69
## 4 1.58 3.27 1.07 0.22 0.56 1.31 2.12 1.88 0.88 0.27
## 5 5.97 11.69 5.89 0.45 0.16 0.01 0.00 0.07 0.06 0.18
## 6 4.41 8.04 4.57 0.77 0.21 0.00 0.02 0.01 0.07 0.05
## Noviembre Diciembre DEF MAM JJA SON Anual est id
## 1 0.64 0.93 1.85 0.94 4.74 3.54 11.07 MODELO.CAMPO.DE.MART 1
## 2 0.43 1.40 5.24 3.46 12.61 3.10 24.41 DONOSO 2
## 3 0.77 0.65 2.40 2.31 5.48 2.96 13.15 VON.HUMBOLDT 3
## 4 0.56 1.39 6.24 1.85 5.31 1.71 15.11 HUAYAN 4
## 5 0.37 2.78 20.44 6.50 0.08 0.61 27.63 SANTA.EULALIA 5
## 6 0.37 1.06 13.51 5.55 0.03 0.49 19.58 CHOSICA 6
## Z Pais
## 1 131 peru
## 2 158 peru
## 3 234 peru
## 4 351 peru
## 5 959 peru
## 6 1220 perurm(list = ls())
#==== Importar a shapefile con 'rgdal'
require(rgdal)
# importar puntos
list.files("Shapefiles", pattern="*.shp")# ver objetos .shp en la carpeta "Shapefiles" ## [1] "cuencas.shp"
## [2] "cuencas.shp.KARLOZ-PC.5896.5652.sr.lock"
## [3] "cuencas.shp.xml"
## [4] "cuencas_est.shp"
## [5] "cuencas_est.shp.xml"
## [6] "ESTACIONES_36.shp"ogrInfo(dsn="Shapefiles", layer="ESTACIONES_36")# propiedades de shape## Source: "Shapefiles", layer: "ESTACIONES_36"
## Driver: ESRI Shapefile; number of rows: 36
## Feature type: wkbPoint with 2 dimensions
## Extent: (256723.4 8652809) - (376289.4 8761413)
## CRS: +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs
## LDID: 89
## Number of fields: 23
## name type length typeName
## 1 FID_ 12 10 Integer64
## 2 x 2 19 Real
## 3 y 2 19 Real
## 4 Enero 2 19 Real
## 5 Febrero 2 19 Real
## 6 Marzo 2 19 Real
## 7 Abril 2 19 Real
## 8 Mayo 2 19 Real
## 9 Junio 2 19 Real
## 10 Julio 2 19 Real
## 11 Agosto 2 19 Real
## 12 Septiembre 2 19 Real
## 13 Octubre 2 19 Real
## 14 Noviembre 2 19 Real
## 15 Diciembre 2 19 Real
## 16 DEF 2 19 Real
## 17 MAM 2 19 Real
## 18 JJA 2 19 Real
## 19 SON 2 19 Real
## 20 Anual 2 19 Real
## 21 est 4 254 String
## 22 id 2 19 Real
## 23 N 2 19 Realdata.shape <- readOGR(dsn="Shapefiles", layer="ESTACIONES_36")# cargar shapefile de estaciones## OGR data source with driver: ESRI Shapefile
## Source: "Shapefiles", layer: "ESTACIONES_36"
## with 36 features
## It has 23 fields
## Integer64 fields read as strings: FID_class(data.shape)# verificar la clase## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"str(data.shape)# verificar si tiene sistema de proyeccion## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 36 obs. of 23 variables:
## .. ..$ FID_ : Factor w/ 1 level "0": 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ x : num [1:36] 279039 256723 287753 269076 318478 ...
## .. ..$ y : num [1:36] 8663781 8731110 8663476 8733416 8681965 ...
## .. ..$ Enero : num [1:36] 0.33 1.07 0.66 1.58 5.97 ...
## .. ..$ Febrero : num [1:36] 0.59 2.77 1.09 3.27 11.69 ...
## .. ..$ Marzo : num [1:36] 0.15 2.15 0.59 1.07 5.89 ...
## .. ..$ Abril : num [1:36] 0.14 0.39 0.67 0.22 0.45 ...
## .. ..$ Mayo : num [1:36] 0.65 0.92 1.05 0.56 0.16 0.21 0.29 0.9 1.32 0.62 ...
## .. ..$ Junio : num [1:36] 1.03 3.61 1.52 1.31 0.01 0 0 0.01 0.2 0 ...
## .. ..$ Julio : num [1:36] 1.25 3.34 1.82 2.12 0 0.02 0 0.01 0 0 ...
## .. ..$ Agosto : num [1:36] 2.46 5.66 2.14 1.88 0.07 0.01 0.01 0.25 0.33 0 ...
## .. ..$ Septiembre: num [1:36] 2.12 2.23 1.5 0.88 0.06 0.07 0.01 1.69 0.6 0.57 ...
## .. ..$ Octubre : num [1:36] 0.78 0.44 0.69 0.27 0.18 0.05 0.29 4.15 8 7.12 ...
## .. ..$ Noviembre : num [1:36] 0.64 0.43 0.77 0.56 0.37 ...
## .. ..$ Diciembre : num [1:36] 0.93 1.4 0.65 1.39 2.78 ...
## .. ..$ DEF : num [1:36] 1.85 5.24 2.4 6.24 20.44 ...
## .. ..$ MAM : num [1:36] 0.94 3.46 2.31 1.85 6.5 ...
## .. ..$ JJA : num [1:36] 4.74 12.61 5.48 5.31 0.08 ...
## .. ..$ SON : num [1:36] 3.54 3.1 2.96 1.71 0.61 ...
## .. ..$ Anual : num [1:36] 11.1 24.4 13.2 15.1 27.6 ...
## .. ..$ est : Factor w/ 36 levels "ANTIOQUIA","ARAHUAY",..: 21 10 35 14 31 9 1 3 18 4 ...
## .. ..$ id : num [1:36] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ N : num [1:36] 1 1 1 1 1 1 1 1 1 1 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:36, 1:2] 279039 256723 287753 269076 318478 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## ..@ bbox : num [1:2, 1:2] 256723 8652809 376289 8761413
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"proj4string(data.shape)# obtener el SCR## [1] "+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"plot(data.shape, axes=T)
# importar polygonos
cuenca.shape <- readOGR(dsn="Shapefiles", layer="cuencas_est")## OGR data source with driver: ESRI Shapefile
## Source: "Shapefiles", layer: "cuencas_est"
## with 3 features
## It has 11 fields
## Integer64 fields read as strings: ID PRINCI_ PRINCI_IDclass(cuenca.shape)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"plot(cuenca.shape, axes=T,col=c("red","cyan","green3"))#col=palette()
plot(data.shape, add=TRUE)
# verificar atributos de SpatialPolygonDataFrame
head(cuenca.shape@data)## ID AREA PERIMETER PRINCI_ PRINCI_ID COD_CUP COD_VTE NOMBRE_P LENGTH
## 0 62 0.19018 2.92220 61 61 1025 1 CHILLON 2.922
## 1 64 0.29771 3.85623 63 63 1026 1 RIMAC 3.856
## 2 66 0.14340 2.52485 65 65 1027 1 LURIN 2.525
## CUENCA Nombre_cue
## 0 1 Cuenca del rÃÂo Chillón
## 1 1 Cuenca del rÃÂo RÃÂmac
## 2 1 Cuenca del rÃÂo LurÃÂn# filtrar los poligonos en base a un valor del atributo
idx <- which(cuenca.shape@data$NOMBRE_P == "RIMAC")
cuenca.rimac <- cuenca.shape[idx,]
class(cuenca.rimac)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"plot(cuenca.rimac, col="cyan", main="Cuenca Rimac")
#==== Proyeccion y Reproyeccion de datos espaciales
# Asignando el sistema de coordenadas
data <- read.csv("Dato_espacial.csv")
coordinates(data) <- ~x + y #alt 126 para "~"
# http://www.spatialreference.org/
# http://www.spatialreference.org/ref/epsg/32718/
proj4string(data)=CRS("+init=epsg:32718")# EPSG:32718: WGS 84 / UTM zone 18S
str(data)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 36 obs. of 20 variables:
## .. ..$ Enero : num [1:36] 0.33 1.07 0.66 1.58 5.97 ...
## .. ..$ Febrero : num [1:36] 0.59 2.77 1.09 3.27 11.69 ...
## .. ..$ Marzo : num [1:36] 0.15 2.15 0.59 1.07 5.89 ...
## .. ..$ Abril : num [1:36] 0.14 0.39 0.67 0.22 0.45 ...
## .. ..$ Mayo : num [1:36] 0.65 0.92 1.05 0.56 0.16 0.21 0.29 0.9 1.32 0.62 ...
## .. ..$ Junio : num [1:36] 1.03 3.61 1.52 1.31 0.01 0 0 0.01 0.2 0 ...
## .. ..$ Julio : num [1:36] 1.25 3.34 1.82 2.12 0 0.02 0 0.01 0 0 ...
## .. ..$ Agosto : num [1:36] 2.46 5.66 2.14 1.88 0.07 0.01 0.01 0.25 0.33 0 ...
## .. ..$ Septiembre: num [1:36] 2.12 2.23 1.5 0.88 0.06 0.07 0.01 1.69 0.6 0.57 ...
## .. ..$ Octubre : num [1:36] 0.78 0.44 0.69 0.27 0.18 0.05 0.29 4.15 8 7.12 ...
## .. ..$ Noviembre : num [1:36] 0.64 0.43 0.77 0.56 0.37 ...
## .. ..$ Diciembre : num [1:36] 0.93 1.4 0.65 1.39 2.78 ...
## .. ..$ DEF : num [1:36] 1.85 5.24 2.4 6.24 20.44 ...
## .. ..$ MAM : num [1:36] 0.94 3.46 2.31 1.85 6.5 ...
## .. ..$ JJA : num [1:36] 4.74 12.61 5.48 5.31 0.08 ...
## .. ..$ SON : num [1:36] 3.54 3.1 2.96 1.71 0.61 ...
## .. ..$ Anual : num [1:36] 11.1 24.4 13.2 15.1 27.6 ...
## .. ..$ est : Factor w/ 36 levels "ANTIOQUIA","ARAHUAY",..: 21 10 35 14 31 9 1 3 18 4 ...
## .. ..$ id : int [1:36] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ Z : int [1:36] 131 158 234 351 959 1220 1881 2164 2442 2523 ...
## ..@ coords.nrs : int [1:2] 1 2
## ..@ coords : num [1:36, 1:2] 279039 256723 287753 269076 318478 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:36] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "x" "y"
## ..@ bbox : num [1:2, 1:2] 256723 8652809 376289 8761413
## .. ..- 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 "+init=epsg:32718 +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"# cambiando de sistema de coordenadas (reproyectando)
data.wgs <- spTransform(data, CRS("+proj=longlat +ellps=WGS84"))
cuenca.wgs <- spTransform(cuenca.shape, CRS("+proj=longlat +ellps=WGS84"))
plot(cuenca.wgs, axes=T, asp=1)
plot(data.wgs, pch=16, col="red", add=TRUE)
# Exportar a google earth
#writeOGR(data.wgs, dsn="Rimac_estaciones.kml", layer="Estacion", driver="KML", dataset_options=c("NameField=est"))#Exportando a google earth
#writeOGR(cuenca.wgs,dsn= "Cuencas.kml",layer= "Cuencas", driver="KML")#Exportando a google earth
#==== Raster
# crear desde cero
mi.raster <- raster(extent(c(-82,-65,-20,0)))
res(mi.raster) <- 2.5#asignar resolucion
mi.raster[] <- 1:ncell(mi.raster)# asignar valores
mi.raster # simple capa## class : RasterLayer
## dimensions : 8, 7, 56 (nrow, ncol, ncell)
## resolution : 2.5, 2.5 (x, y)
## extent : -82, -64.5, -20, 0 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 1, 56 (min, max)summary(mi.raster)## layer
## Min. 1.00
## 1st Qu. 14.75
## Median 28.50
## 3rd Qu. 42.25
## Max. 56.00
## NA's 0.00plot(mi.raster)
# algebra u operaciones matematicas en raster
mi.raster2 <- -1*mi.raster^2
plot(mi.raster2)
summary(mi.raster2)## layer
## Min. -3136.00
## 1st Qu. -1785.25
## Median -812.50
## 3rd Qu. -217.75
## Max. -1.00
## NA's 0.00# Importar un raster
lista <- list.files("Pisco_clima_mensual", pattern=".tif", full.names=TRUE)
lista## [1] "Pisco_clima_mensual/01.tif" "Pisco_clima_mensual/02.tif"
## [3] "Pisco_clima_mensual/03.tif" "Pisco_clima_mensual/04.tif"
## [5] "Pisco_clima_mensual/05.tif" "Pisco_clima_mensual/06.tif"
## [7] "Pisco_clima_mensual/07.tif" "Pisco_clima_mensual/08.tif"
## [9] "Pisco_clima_mensual/09.tif" "Pisco_clima_mensual/10.tif"
## [11] "Pisco_clima_mensual/11.tif" "Pisco_clima_mensual/12.tif"GDALinfo(lista[1])#ver informacion acerca del raster## rows 400
## columns 267
## bands 1
## lower left origin.x -81.35
## lower left origin.y -19
## res.x 0.05
## res.y 0.05
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=longlat +datum=WGS84 +no_defs
## file Pisco_clima_mensual/01.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -1.7e+308 3 267
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 0 730.7213 160.6065 127.635
## Metadata:
## AREA_OR_POINT=Areapp.ene.1 <- readGDAL(lista[1])# utilizando el paquete GDAL## Pisco_clima_mensual/01.tif has GDAL driver GTiff
## and has 400 rows and 267 columnsstr(pp.ene.1)# Que clase de objeto es## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
## ..@ data :'data.frame': 106800 obs. of 1 variable:
## .. ..$ band1: num [1:106800] 99 99.3 97 100.5 99.5 ...
## ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
## .. .. ..@ cellcentre.offset: Named num [1:2] -81.3 -19
## .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
## .. .. ..@ cellsize : num [1:2] 0.05 0.05
## .. .. ..@ cells.dim : int [1:2] 267 400
## ..@ bbox : num [1:2, 1:2] -81.3 -19 -68 1
## .. ..- 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 +ellps=WGS84 +towgs84=0,0,0"pp.ene.2 <- raster(lista[1])# utilizando el paquete raster
str(pp.ene.2)# Que clase de objeto es## Formal class 'RasterLayer' [package "raster"] with 12 slots
## ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots
## .. .. ..@ name : chr "E:\\SENAMHI_2017\\Curso_R_SENAMHI\\3_Datos_espaciales\\Pisco_clima_mensual\\01.tif"
## .. .. ..@ datanotation: chr "FLT8S"
## .. .. ..@ byteorder : chr "little"
## .. .. ..@ nodatavalue : num -Inf
## .. .. ..@ NAchanged : logi FALSE
## .. .. ..@ nbands : int 1
## .. .. ..@ bandorder : chr "BIL"
## .. .. ..@ offset : int 0
## .. .. ..@ toptobottom : logi TRUE
## .. .. ..@ blockrows : int 3
## .. .. ..@ blockcols : int 267
## .. .. ..@ driver : chr "gdal"
## .. .. ..@ open : logi FALSE
## ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots
## .. .. ..@ values : logi(0)
## .. .. ..@ offset : num 0
## .. .. ..@ gain : num 1
## .. .. ..@ inmemory : logi FALSE
## .. .. ..@ fromdisk : logi TRUE
## .. .. ..@ isfactor : logi FALSE
## .. .. ..@ attributes: list()
## .. .. ..@ haveminmax: logi TRUE
## .. .. ..@ min : num 0
## .. .. ..@ max : num 731
## .. .. ..@ band : int 1
## .. .. ..@ unit : chr ""
## .. .. ..@ names : chr "X01"
## ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots
## .. .. ..@ type : chr(0)
## .. .. ..@ values : logi(0)
## .. .. ..@ color : logi(0)
## .. .. ..@ names : logi(0)
## .. .. ..@ colortable: logi(0)
## ..@ title : chr(0)
## ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots
## .. .. ..@ xmin: num -81.3
## .. .. ..@ xmax: num -68
## .. .. ..@ ymin: num -19
## .. .. ..@ ymax: num 1
## ..@ rotated : logi FALSE
## ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
## .. .. ..@ geotrans: num(0)
## .. .. ..@ transfun:function ()
## ..@ ncols : int 267
## ..@ nrows : int 400
## ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
## ..@ history : list()
## ..@ z : list()# plot
plot(pp.ene.2)
# datos basicos del raster
ncell(pp.ene.2)# numero de celdas## [1] 106800extent(pp.ene.2)# limites## class : Extent
## xmin : -81.35
## xmax : -68
## ymin : -19
## ymax : 0.9999993names(pp.ene.2)# nombre raster## [1] "X01"ncol(pp.ene.2)# numero de columnas## [1] 267nrow(pp.ene.2)# numero filas## [1] 400# valores del raster
a <- pp.ene.2@data@values# que paso
b <- pp.ene.2[]
c <- getValues(pp.ene.2)
identical(b,c)## [1] TRUErm(a,b,c)
# Importar multiples archivos de raster
r <- stack(lista)
class(r)## [1] "RasterStack"
## attr(,"package")
## [1] "raster"summary(r)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (93.63% of all cells)## X01 X02 X03 X04 X05 X06
## Min. 0.000000 0.00000 0.00000 0.2027299 0.000000 0.000000
## 1st Qu. 7.357911 10.52383 7.03367 7.0138608 3.274348 4.493996
## Median 166.424207 179.20131 237.70641 159.7383888 83.754601 46.636537
## 3rd Qu. 262.633921 279.11331 340.91202 294.8741512 214.594558 142.562478
## Max. 730.721259 926.94256 827.22566 665.5253296 579.017504 501.457341
## NA's 0.000000 0.00000 0.00000 0.0000000 0.000000 0.000000
## X07 X08 X09 X10 X11 X12
## Min. 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## 1st Qu. 2.120191 4.133392 3.104162 3.475238 4.316219 2.054533
## Median 38.620461 40.315495 70.737636 121.280697 132.671229 151.995335
## 3rd Qu. 115.013021 109.964138 153.097635 205.830888 235.944032 272.011611
## Max. 558.355537 390.139224 390.452438 602.050851 820.022662 919.821390
## NA's 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000plot(r)
projection(r)## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"# br <- brick(r)
# class(br)
# Cortar raster para las cuencas de CHIRILU
plot(r[[1]])
plot(cuenca.wgs, add=T)
r.chirilu <- crop(r, cuenca.wgs, snap="out")
plot(r.chirilu) 
r.chirilu <- mask(r.chirilu, cuenca.wgs)
#plot(r.chirilu)
names(r.chirilu) ## [1] "X01" "X02" "X03" "X04" "X05" "X06" "X07" "X08" "X09" "X10" "X11"
## [12] "X12"names(r.chirilu) <- month.abb
plot(r.chirilu) 
# plot utilizando paquete lattice y latticeExtra
library(lattice)
library(latticeExtra)
spplot(r.chirilu, col.regions = terrain.colors(100))
spplot(r.chirilu, col.regions = rev(terrain.colors(100)))
range(r.chirilu)## class : RasterBrick
## dimensions : 20, 24, 480, 2 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -77.2, -76, -12.3, -11.3 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : range_min, range_max
## min values : 0.00569161, 1.40154847
## max values : 13.5527, 152.1218spplot(r.chirilu,
col.regions = rev(terrain.colors(100)),
at=seq(0,160,5),
scales = list(draw=TRUE)
) +
layer(sp.polygons(cuenca.wgs, lwd=0.8, col='darkgray'))+
layer(sp.points(data.wgs, lwd=0.8, col='red'))
# extraer datos desde raster
# para puntos
#pp.estacion <- extract(r.chirilu,data.wgs)
pp.estacion <- extract(r,data.wgs)
row.names(pp.estacion) <- data.wgs@data$est
colnames(pp.estacion) <- month.abb
View(pp.estacion)
boxplot(pp.estacion, main="Variabilidad y estacionalidad de precipitacion - CHIRILU [mm]")
# Para poligonos
#pp.cuenca <- extract(r.chirilu,cuenca.wgs,fun=mean)
pp.cuenca <- extract(r,cuenca.wgs,fun=mean)
row.names(pp.cuenca) <- cuenca.wgs@data$NOMBRE_P
colnames(pp.cuenca) <- month.abb
View(pp.cuenca)
plot(pp.cuenca[1,], type="o", col="1", ylim=c(0,110), ylab="Prec. [mm]", xlab = "Meses", main="Prec. prom areal - CHIRILU [mm]")
lines(pp.cuenca[2,], type="o", col="2")
lines(pp.cuenca[3,], type="o", col="3")
legend("topright", as.character(cuenca.wgs@data$NOMBRE_P), lty = 1, col = c(1:3),cex=0.75)
# ------------------------------------------------------------------------
# LEER DATOS DE PISCO - PRECIPITACION MENSUAL DEL PERIODO 1981 - 2016-----
# ------------------------------------------------------------------------
#install.packages("ncdf4")
library(ncdf4)
library(raster)
#Pisco.prec.stack <- stack("PISCOV3-MONTHLY.nc")# leer netcdf con stack
#Pisco.prec.stack
Pisco.prec.brick <- brick("PISCOV3-MONTHLY.nc")# leer netcdf con brick
Pisco.prec.brick## class : RasterBrick
## dimensions : 400, 267, 106800, 432 (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05 (x, y)
## extent : -81.35, -68, -19, 0.9999993 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : E:\SENAMHI_2017\Curso_R_SENAMHI\3_Datos_espaciales\PISCOV3-MONTHLY.nc
## names : X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, ...
## months since 1981-01-01 12:00:00: 0, 431 (min, max)
## varname : variablenlayers(Pisco.prec.brick)## [1] 432plot(Pisco.prec.brick[[1:12]])
# Extraer los valores prec. promedio areal para las cuencas CHIRILU
pp.cuenca.mensual <- extract(Pisco.prec.brick, cuenca.wgs, fun=mean)
row.names(pp.cuenca.mensual) <- cuenca.wgs@data$NOMBRE_P
colnames(pp.cuenca.mensual) <- 1:ncol(pp.cuenca.mensual)
View(pp.cuenca.mensual)
range(pp.cuenca.mensual)## [1] 0.04179044 157.59823268plot(pp.cuenca.mensual[1,], type="o", col="1", ylim=c(0,200), ylab="Prec. [mm]", xlab = "Meses", main="Prec. prom areal - CHIRILU [mm]")
lines(pp.cuenca.mensual[2,], type="o", col="2")
lines(pp.cuenca.mensual[3,], type="o", col="3")
legend("topright", as.character(cuenca.wgs@data$NOMBRE_P), lty = 1, col = c(1:3),cex=0.75)
# ------------------------------------------------------------------------
# Importar datos vectoriales Online con 'raster'---------------------------
# ------------------------------------------------------------------------
library(raster)
args(getData)## function (name = "GADM", download = TRUE, path = "", ...)
## NULL#getData('ISO3')# Codigo ISO de los paises
# Shapefile de limites departamentales del peru
peru.dptos <- getData('GADM', country='PER', level=1)
class(peru.dptos)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"plot(peru.dptos, axes=T, asp=1)
# datos de elevacion para todo el peru
dem <- getData('alt', country='PER')[[1]]#dem peru
#str(dem)
plot(dem)
# datos de climatologia de precipitacion mensual "worldclim"
prec <- getData('worldclim', var='prec', res=0.5, lon=-77, lat=-12)
plot(prec)
# plot prec enero + shape de departamentos
plot(prec[[1]])
plot(peru.dptos, add=TRUE)
# FIN