Curso de Geoestadística.
Magíster en Teledetección.
Universidad Mayor.
Estudiante: Nelson Mattie
Profesor: Francisco Zambrano
Fecha: 13 de febrero de 2021
1 Introducción.
En el presente informe investigamos las principales técnicas de geoestadística y las aplicamos a un conjunto de datos geográficos en una zona de interés con el propósito de asistir en la toma de decisiones en ámbitos de planificación territorial, investigación, conservación, meteorología, entre otros. Contamos con datos geográficos de tipo meteorológico disponibles inicialmente a una extensión planetaria, los que fuimos sucesivamente procesando para identificar zonas de menor tamaño que fueran mayor interés para nosotros.
Usamos técnicas tales como variograma que nos permite investigar la estructura de dependencia en las variables geográficas que disponemos. Conocer si las variables presentan dependencia es importante porque es una propiedad que da respaldo al momento de construir un modelo de predicción espacial.
2 Área de estudio.
2.0.1 Configurar documento.
knitr::knit_hooks$set(time_it = local({
now <- NULL
function(before, options) {
if (before) {
# record the current time before each chunk
now <<- Sys.time()
} else {
# calculate the time difference after a chunk
res <- difftime(Sys.time(), now)
# return a character string to show the time
print(paste("Time for this code chunk to run:", res))
}
}
}))
knitr::opts_chunk$set(time_it = TRUE)2.0.2 Preparar el entorno de trabajo.
ruta_proyecto <- dirname(getwd())
ruta_clase_01 <- file.path(ruta_proyecto, "Taller1")
ruta_clase_02 <- file.path(ruta_proyecto, "Taller2")
ruta_clase_03 <- file.path(ruta_proyecto, "Taller3")
ruta_clase_04 <- file.path(ruta_proyecto, "Taller4")
ruta_datos <- file.path(ruta_proyecto, "data")if(!require(tmap)) suppressMessages( suppressWarnings( install.packages("tmap") ))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))
chile_zona_a <- sf::st_read(file.path(ruta_datos,"shp/zona_a.shp"),quiet=TRUE)
chile_zona_a <- sf::st_as_sf(chile_zona_a)
tm_shape(chile_zona_a)+tm_polygons("NAME_1", legend.show = F) + tm_text("NAME_1")
3 Materiales y Métodos.
3.1 Datos.
3.1.1 Preparar los datos.
3.1.1.1 Explorar zonas dentro de Chile para el análisis.
Podemos obtener un mapa de Chile con las divisiones político administrativas regionales en donde podemos situar con más detalle las coordenadas de las estaciones de temperaturas máximas.
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
if(file.exists(file.path(ruta_datos,"rds/gadm36_CHL_1_sp.rds")))
{ limite_chile <- readRDS(file.path(ruta_datos,"rds/gadm36_CHL_1_sp.rds"))
} else {
limite_chile <- getData('GADM',country='chl',level=1)
}
class(limite_chile)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"Del mismo modo, lo convertimos a un objeto de tipo sf.
limite_chile_sf <- st_as_sf(limite_chile)
class(limite_chile_sf)## [1] "sf" "data.frame"Aprovechamos de exportar los límites regionales de Chile en formato shapefile para usarlo en caso de no disponer de internet o para operaciones con herramientas externas a modo de verificar consistencia con cada geoproceso dentro de R.
if(!require(rgdal)) suppressMessages( suppressWarnings( install.packages("rgdal") ))
writeOGR(as(limite_chile_sf,'Spatial'),
layer="limite_chile_sf",
driver="ESRI Shapefile",
layer_options = "ENCODING=UTF-8",
file.path(ruta_datos,"shp/limite_chile_sf.shp"), overwrite=TRUE)Otra manera puede ser leyendo desde un archivo shapefile en disco. Observamos una gráfica del mapa de Chile.
png(file.path(ruta_datos,"figs/limite_chile_sf.png"), units="in", width=6, height=6, res=300)
plot(st_geometry(limite_chile_sf))
dev.off()## png
## 2
Cargamos cartografía de Chile para definir algunas zonas de interés.
if(!require(rgdal)) suppressMessages( suppressWarnings( install.packages("rgdal") ))
limite_chile_sf <- st_read(file.path(ruta_datos,"shp/limite_chile_sf.shp"),quiet=TRUE)
zona_a <- limite_chile_sf [ c(8,12),]
writeOGR(as(zona_a,'Spatial'),
layer="zona_a",
driver="ESRI Shapefile",
layer_options = "ENCODING=UTF-8",
file.path(ruta_datos,"shp/zona_a.shp"), overwrite=TRUE)Seleccionamos dos regiones para graficarlas y que serán nuestra área de estudio.
limite_chile_a <- limite_chile_sf[c(8,12),]
png(file.path(ruta_datos,"figs/limite_chile_zona_a.png"), units="in", width=6, height=6, res=300)
plot(st_geometry(limite_chile_a))
dev.off()## png
## 2
Creamos la zona de estudio para esto unimos las dos regiones recién seleccionadas en una sola región.
zona_a <- st_union(limite_chile_a) #area de estudio## although coordinates are longitude/latitude, st_union assumes that they are planar3.1.1.2 Explorar la variable Tmax.
Los datos se encuentran en un objeto data.frame que fue almacenado en un archivo RDS que sirve para guardar objetos de forma individual en disco. Se usa la función readRDS y si están en .csv read.csv2.
tmax_2017 <- readRDS(file.path(ruta_datos,"rds/tmax_2017_taller1.rds"))La variable Tmax está medida para los meses Enero, Abril, Julio y Octubre, por separado, es decir, tenemos cuatro variables Tmax que presenta cada mes distinto. Miramos las observaciones faltantes para los datos tmax_2017.
faltantes_tmax_2017 <- apply( tmax_2017, 2, function(x) sum(is.na(x)))
faltantes_tmax_2017## codigo_estacion institucion altura latitud longitud Jan2017 Apr2017 Jul2017 Oct2017
## 0 0 0 0 0 544 551 565 559Tenemos que existen observaciones faltantes en ciertas estaciones para cada fecha.
pivot_codigo_estacion <- aggregate( institucion ~ codigo_estacion, data=tmax_2017, FUN=length )
pivot_codigo_estacion[pivot_codigo_estacion['institucion']>1,]## [1] codigo_estacion institucion
## <0 rows> (or 0-length row.names)nrow(pivot_codigo_estacion)## [1] 909Observamos que no existen estaciones con registros repetidos en la base de datos actual con un total de 909 estaciones distintas, lo que equivale al número total de observaciones.
class(tmax_2017)## [1] "tbl_df" "tbl" "data.frame"Podemos crear un objeto sf a partir de un data.frame indicando explícitamente cuales son las coordenadas.
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))
if(!require(spData)) suppressMessages( suppressWarnings( install.packages("spData") ))tmax_2017_sf <- st_as_sf(tmax_2017,coords = c('longitud','latitud'))
st_crs(tmax_2017_sf) <- 4326
class(tmax_2017_sf)## [1] "sf" "tbl_df" "tbl" "data.frame"st_crs(tmax_2017_sf)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4326"]]Con lo cual tenemos un objeto del tipo simplefeature a partir de un data.frame.
png(file.path(ruta_datos,"figs/mundo_estaciones.png"), units="in", width=6, height=6, res=300)
plot(st_geometry(world))
plot(tmax_2017_sf["altura"],add=TRUE)
dev.off()## png
## 2
3.1.1.3 Explorar la variable TSS de los datos raster MODIS.
3.1.1.3.1 Usando las imágenes obtenidas inicialmente.
Optamos por renombrar los archivos raster para abreviar los nombres al momento de graficar, manteniendo siempre copias de los archivos originales.
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_001.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_001.tif"), file.path(ruta_datos,"raster/MOD11C3_2017_001.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_091.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_091.tif"), file.path(ruta_datos,"raster/MOD11C3_2017_091.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_182.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_182.tif"), file.path(ruta_datos,"raster/MOD11C3_2017_182.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_274.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2017_274.tif"), file.path(ruta_datos,"raster/MOD11C3_2017_274.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_001.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_001.tif"), file.path(ruta_datos,"raster/MOD11C3_2020_001.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_122.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_122.tif"), file.path(ruta_datos,"raster/MOD11C3_2020_122.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_183.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_183.tif"), file.path(ruta_datos,"raster/MOD11C3_2020_183.tif"))
if(file.exists(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_306.tif"))) file.rename(file.path(ruta_datos,"raster/MOD11C3_LST_Day_CMG_2020_306.tif"), file.path(ruta_datos,"raster/MOD11C3_2020_306.tif"))Cargamos los archivos raster 2017 en modo stack desde la carpeta raster.
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
# mod11 <- raster(file.path(ruta_datos,"raster/MOD11C3_2017_001.tif")) #para una sola capa
archivos_mod11 <- list.files(file.path(ruta_datos,"raster"),pattern='MOD11C3_2017_\\d+.tif$')
archivos_mod11## [1] "MOD11C3_2017_001.tif" "MOD11C3_2017_091.tif" "MOD11C3_2017_182.tif" "MOD11C3_2017_274.tif"stack_mod11 <- stack(file.path(ruta_datos,"raster",archivos_mod11))
stack_mod11## class : RasterStack
## dimensions : 8050, 16100, 129605000, 4 (nrow, ncol, ncell, nlayers)
## resolution : 2486.243, 2486.088 (x, y)
## extent : -20015109, 20013395, -10005454, 10007555 (xmin, xmax, ymin, ymax)
## crs : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
## names : MOD11C3_2017_001, MOD11C3_2017_091, MOD11C3_2017_182, MOD11C3_2017_274
## min values : 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535Asignamos SRC modis.
crs(stack_mod11) <- CRS("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext")Las imágenes originales no contenían un sistema de proyección, por eso es que lo definimos dentro de la sesión R. Exportaremos una copia de cada imagen en disco con la proyección recién definida.
Cambiamos SRC de la zona de estudio.
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))
zona_a_modis <- st_transform(zona_a,crs(stack_mod11))
stack_mod11 <- crop(stack_mod11,as(zona_a_modis,'Spatial'))Este es el punto en donde optamos por cambiar de estrategia para procesar nuestros datos geográficos, puesto que no creemos buena idea realizar los geoprocesos con el sistema de proyección sinusoidal modis para todas las entidades geográficas en vez del sistema WGS84 que es el más ampliamente usado.
writeRaster(stack_mod11, file.path(ruta_datos,"raster/mod11sinu.EAJO2017.tif"), overwrite=TRUE)Volveremos a generar esta variable desde nuevas imágenes descargadas desde el proveedor NASA como explicaremos a continuación.
3.1.1.3.2 Usando imágenes descargadas del proveedor.
Debido a que asignamos sistemas de proyección y reproyectamos varias veces las imágenes MODIS inicialmente recibidas llegamos a la conclusión de que éstas no reflejan correctamente, ni las unidades en que queremos expresarlas, ni el sistema de proyección que queremos trabajar con el resto de las entidades geográficas. Por este motivo optamos por ingresar al sitio web del proveedor, crear una cuenta y descargar nuevamente las imágenes para procesarlas desde cero en caso de que fuera este el motivo de que éstas no reflejaran las propiedades que quisimos asignarles.
Nos ayudamos de QGis y la librería GDAL para ir cotejando los rasters con lo procesado en R. Investigamos en el sitio web NASA las capas que contiene cada imagen para seleccionar la que corresponde a Land Surface Temperature o LST.
Luego de crear una cuenta exploramos el repositorio de imágenes MODIS primero por año y luego seleccionamos los meses EAJO. Los meses están indicados por el número de día correlativo en el año. Las imágenes están en formato hdf4 y cada una contiene 17 capas de las cuales es la primera la que contienen la variable LST o TSS para nuestro estudio.
Una vez descargadas las imágenes de 2017, extrajimos la primera capa desde cada archivo hdf4 para los cuatro meses y las proyectamos con EPSG 4326 para tener el sistema WGS84 y sean compatibles con las demás entidades geográficas.
Generaremos las mismas imágenes con proyección WGS84 (EPSG 4326) para ver luego con cuál de ellas resulta mejor trabajar. Hemos visto que hay dos razones por las cuales convendría seleccionar de nuevo el producto MODIS para este trabajo.
3.2 Diagrama.
el producto MODIS para este trabajo. La primera es porque cada imagen es del mapa mundi completo, por lo que hay mucha información que no necesitamos. La segunda es porque la proyección sinusoidal es muy poco adecuada para trabajar con las otras entidades geográficas del estudio. Por estos motivos veremos algunas opciones avanzadas para consultar el proveedor de imágenes modis usando el paquete R MODIS.
Cargamos nuestras nuevas imágenes modis para trabajar con la variable TSS. Verificamos que contenga el sistema WGS84.
archivos_modis_tss <- list.files(file.path(ruta_datos,"raster"),pattern="MOD11C3_A2017\\d+.tif")
archivos_modis_tss## [1] "MOD11C3_A2017001.tif" "MOD11C3_A2017091.tif" "MOD11C3_A2017182.tif" "MOD11C3_A2017274.tif"modis_tss_stack <- stack(file.path(ruta_datos,"raster",archivos_modis_tss))
crs(modis_tss_stack)## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,03.2.0.1 Diagrama de proceso de los datos.
En esta sección y en las siguientes realizamos varios procesos con los datos, incluyendo geoprocesos tanto con el programa R como con algunas herramientas de apoyo. Conceptualmente es importante contar con un diagrama que nos permita diseñar los pasos que realizaremos sobre la información ya que es sabido que un ambiente cambiante, es clave contar con una hoja de ruta que nos permita mantener la claridad en las operaciones requeridas con entidades espaciales.
El diagrama anterior muestra que primero se trabajará con las variables para realizar un análisis exploratorio, principalmente en forma consolidada. Es decir, las variables serán consolidadas de manera
4 Resultados.
4.1 Análisis Exploratorio.
consolidada. Es decir, las variables serán consolidadas de manera que adquieran la estructura de la variable dependiente Tmax, que corresponde a puntos espaciales. A pesar que se también se las analizará en forma individual, la mayoría de los análisis se realizan sobre una misma estructura de datos.
Posteriormente, los pasos de regresión lineal, variograma, se trabajará con las variables consolidadas usando la estructura de puntos de la variable Tmax. Finalmente, la etapa de predicción espacial, aunque la predicción se puede hacer, tanto en forma de coordenadas como de una cobertura continua, representada en este caso, por una imagen satelital, es esta última estructura la que usaremos para realizar nuestras predicciones.
4.1.0.1 Intersectar variable TSS con Tmax para Chile.
Usar un poco el paquete {tidyr} para ordenar los datos.
if(!require(tidyr)) suppressMessages( suppressWarnings( install.packages("tidyr") ))
if(!require(dplyr)) suppressMessages( suppressWarnings( install.packages("dplyr") ))Dejaremos disponibles datos intersectados entre TSS con Tmax para todo Chile en caso de que posteriormente hagamos uso de ellos. Primero convetimos Tmax a formato longer para descartar posteriormente los datos faltantes e intersectar con TSS.
tmax_2017_sf_nocoord <- st_drop_geometry(tmax_2017_sf)
tmax_2017_sf_longer <- pivot_longer(tmax_2017_sf_nocoord, 4:7,names_to='Meses',values_to = 'Tmax')Filtramos las observaciones faltantes para la variable Tmax que se encuentra en formato longer.
tmax_2017_sf_longer_filter <- tmax_2017_sf_longer %>%
filter(!is.na(Tmax))
dim(tmax_2017_sf_longer_filter)## [1] 1417 5En el siguiente paso filtramos las estaciones que solamente se encuentran en Chile.
tmax_2017_sf_estaciones <- filter(tmax_2017_sf,
codigo_estacion %in%
tmax_2017_sf_longer_filter$codigo_estacion)
dim(tmax_2017_sf_estaciones)## [1] 382 8Intersectamos Tmax con TSS para todo Chile. Lo que obtenemos son los datos TSS solamente a nivel de Chile, pero no los estamos uniendo a los datos Tmax en este paso.
tmax_tss_2017_chile_extraccion <- raster::extract(modis_tss_stack,tmax_2017_sf_estaciones)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
tmax_tss_2017_chile_extraccion <- apply(tmax_tss_2017_chile_extraccion,2,FUN = function(x) x - 273.15)Unimos los datos TSS con Tmax para posterior análisis a nivel de Chile en formato wide.
tmax_tss_2017_chile_unidos <- cbind(tmax_2017_sf_estaciones,tmax_tss_2017_chile_extraccion)Contamos cuantas observaciones por cada mes tenemos disponibles.
aggregate(Tmax~Meses,tmax_2017_sf_longer,length)## Meses Tmax
## 1 Apr2017 358
## 2 Jan2017 365
## 3 Jul2017 344
## 4 Oct2017 3504.1.0.2 Análisis exploratorio en la zona a.
Filtramos los datos de estaciones en la zona de estudio.
tmax_2017_sf_zona_a <- tmax_2017_sf[zona_a,]## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planarAnálisis exploratorio de los datos tabulados el paquete {ggplot2}.
dim(tmax_2017_sf_zona_a)## [1] 37 8Separamos los datos de sus coordenadas para trabajar solamente con los datos.
tmax_2017_sf_zona_a_nocoord <- st_drop_geometry(tmax_2017_sf_zona_a)
tmax_2017_sf_zona_a_longer <- pivot_longer(tmax_2017_sf_zona_a_nocoord, 4:7,names_to='Meses',values_to = 'Tmax')Renombrar meses para conservar orden original de cada mes y evitar el orden alfabético en los sucesivos gráficos.
tmax_2017_sf_zona_a_longer$Meses [tmax_2017_sf_zona_a_longer$Meses == 'Jan2017'] <- '2017_01'
tmax_2017_sf_zona_a_longer$Meses [tmax_2017_sf_zona_a_longer$Meses == 'Apr2017'] <- '2017_04'
tmax_2017_sf_zona_a_longer$Meses [tmax_2017_sf_zona_a_longer$Meses == 'Jul2017'] <- '2017_07'
tmax_2017_sf_zona_a_longer$Meses [tmax_2017_sf_zona_a_longer$Meses == 'Oct2017'] <- '2017_10'
aggregate( Tmax ~ Meses , data=tmax_2017_sf_zona_a_longer, FUN=length)## Meses Tmax
## 1 2017_01 23
## 2 2017_04 23
## 3 2017_07 23
## 4 2017_10 23Graficar los datos usando el paquete ggplot2.
if(!require(ggplot2)) suppressMessages( suppressWarnings( install.packages("ggplot2") ))ggplot(tmax_2017_sf_zona_a_longer) +
geom_boxplot(aes(Meses,Tmax)) +
theme_minimal()## Warning: Removed 56 rows containing non-finite values (stat_boxplot).
ggplot(tmax_2017_sf_zona_a_longer) +
geom_density(aes(Tmax)) +
facet_wrap(~Meses) +
theme_minimal()## Warning: Removed 56 rows containing non-finite values (stat_density).
ggsave(file.path(ruta_datos,"figs/density_esta.png"))## Saving 7 x 5 in image## Warning: Removed 56 rows containing non-finite values (stat_density).Los datos faltantes en la zona de estudio corresponden a 56 observaciones, mientras que los datos presentes son 92, con 23 observaciones para cada mes en cuestión.
dim(tmax_2017_sf_zona_a_longer)## [1] 148 5dim(tmax_2017_sf_zona_a_longer[is.na(tmax_2017_sf_zona_a_longer$Tmax),])## [1] 56 5dim(tmax_2017_sf_zona_a_longer[!is.na(tmax_2017_sf_zona_a_longer$Tmax),])## [1] 92 5tmax_2017_sf_zona_a_longer %>%
group_by(Meses) %>%
summarize(n=sum(!is.na(Tmax)))## # A tibble: 4 x 2
## Meses n
## * <chr> <int>
## 1 2017_01 23
## 2 2017_04 23
## 3 2017_07 23
## 4 2017_10 23Excluimos los datos faltantes en un nuevo data.frame.
tmax_2017_sf_zona_a_longer_filter <- tmax_2017_sf_zona_a_longer %>%
filter(!is.na(Tmax))
dim(tmax_2017_sf_zona_a_longer_filter)## [1] 92 5Graficamos los datos tmax filtrados para la zona de interés.
ggplot(tmax_2017_sf_zona_a_longer_filter) +
geom_boxplot(aes(Meses,Tmax)) +
theme_minimal()
ggplot(tmax_2017_sf_zona_a_longer_filter) +
geom_density(aes(Tmax)) +
labs(caption = 'n = 23') +
facet_wrap(~Meses) +
theme_minimal()
Ahora usando la zona de estudio con la proyección 4326 inicial, cortamos la variable TSS para acotar el área a la zona correspondiente a la región de O'Higgins y del Maule, lo que denominamos zona_a.
modis_tss_stack_zona_a <- crop(modis_tss_stack, as(zona_a,'Spatial'))Grabaremos una copia de la zona a para uso posterior.
saveRDS( zona_a, file.path(ruta_datos,"rds/chile_zona_a.rds"))4.1.0.3 Intersectar variable TSS con Tmax para la zona a.
Ahora tenemos los datos raster y vamos a extraer los puntos de las estaciones filtradas.
Filtramos los datos de tmax_2017_sf_zona_a que tienen geometría.
datos_estaciones <- filter(tmax_2017_sf_zona_a,
codigo_estacion %in%
tmax_2017_sf_zona_a_longer_filter$codigo_estacion)
dim(datos_estaciones)## [1] 23 8datos_extraidos <- raster::extract(modis_tss_stack_zona_a,datos_estaciones)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
datos_extraidos <- apply(datos_extraidos,2,FUN = function(x) x - 273.15)
dim(datos_extraidos)## [1] 23 4Datos estaciones + datos modis con datos filtrados y geometría.
data_unida <- cbind(datos_estaciones,datos_extraidos)
dim(data_unida)## [1] 23 12datos_sin_geometria <- st_drop_geometry(data_unida)plot(datos_sin_geometria$Jan2017,datos_sin_geometria$MOD11C3_2017_001)
Transformamos los datos para crear un ggplot.
datos_sin_geometria1 <- pivot_longer(datos_sin_geometria,4:7,names_to = 'Meses',values_to = 'Tmax_esta')
datos_sin_geometria2 <- pivot_longer(datos_sin_geometria,8:11,names_to = 'Meses', values_to = 'Tmax_modis')Datos unidos en formato pivot_longer.
dataFinal <- cbind(datos_sin_geometria1,datos_sin_geometria2)
dataFinal <- dataFinal[,c(1:3,8,9,18)]4.1.0.4 Visualizar datos de TSS.
ggplot(dataFinal) +
geom_point(aes(Tmax_modis,Tmax_esta)) +
facet_wrap(~Meses,scales = 'free') +
theme_bw()
4.1.0.5 Visualizar datos de Tmax.
Haremos algunas visualizaciones de datos de estaciones usando paquete {tmap}.
if(!require(tmap)) suppressMessages( suppressWarnings( install.packages("tmap") ))mapEst1 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='Jan2017',title='Enero 2017')mapEst2 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='Apr2017',title='Abril 2017')mapEst3 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='Jul2017',title='Julio 2017')mapEst4 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='Oct2017',title='Octubre 2017')mapsEsta <- tmap_arrange(mapEst1,mapEst2,mapEst3,mapEst4,ncol=2)
tmap::tmap_options(show.messages=FALSE)
tmap_save(mapsEsta,file.path(ruta_datos,"figs/var_espacial_tmax_estaciones.png"))4.1.0.6 Mapa de variación espacial de temperatura de MODIS.
names(data_unida)[8:11] <- c('TmodisJan','TmodisApr','TmodisJul','TmodisOct')mapMod1 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='TmodisJan',title='Enero 2017')mapMod2 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='TmodisApr',title='Abril 2017')mapMod3 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='TmodisJul',title='Julio 2017')mapMod4 <- tm_shape(zona_a) +
tm_borders() +
tm_shape(data_unida) +
tm_dots(shape=5,size = .3,col='TmodisOct',title='Octubre 2017')Guardamos el data.frame con los datos unidos en el disco.
saveRDS(data_unida,file.path(ruta_datos,"rds/data_unida_tmax_mod11_2017.rds"))mapsEsta <- tmap_arrange(mapMod1,mapMod2,mapMod3,mapMod4,ncol=2)
tmap::tmap_options(show.messages=FALSE)
tmap_save(mapsEsta,file.path(ruta_datos,"figs/var_espacial_tmax_MOD11.png"))Todos los mapas juntos estaciones y MOD11.
mapsTodos <- tmap_arrange(mapEst1,mapEst2,mapEst3,mapEst4,
mapMod1,mapMod2,mapMod3,mapMod4,nrow=2)modis_tss_stack_zona_a_masc <- mask(modis_tss_stack_zona_a,as(zona_a,'Spatial'))Graficamos resultado en zona de interés.
png(file.path(ruta_datos,"figs/modis_tss_stack_zona_a_masc.png"), units="in", width=6, height=6, res=300)
plot(modis_tss_stack_zona_a_masc)
dev.off()## png
## 2
Grabamos el stack con las imágenes raster de la zona del estudio en disco.
# crs(modis_tss_stack) <- sp::CRS('+init=EPSG:4326')
# crs(modis_tss_stack) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
writeRaster(modis_tss_stack, file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack.tif"), overwrite=TRUE)Grabamos el raster con la zona de estudio para uso posterior.
writeRaster(modis_tss_stack_zona_a, file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack_zona_a.tif"), overwrite=TRUE)4.1.1 Resultados del taller 1.
source(file.path(ruta_proyecto,"twee.R"))4.1.1.1 Gráficos y mapas grabados en formato png.
twee(file.path(ruta_proyecto,"anexos/taller1/png"))## -- density_esta.png
## -- diagrama.png
## -- limite_chile_sf.png
## -- limite_chile_zona_a.png
## -- modis_tss_stack_zona_a_masc.png
## -- mundo_estaciones.png
## -- var_espacial_tmax_estaciones.png
## -- var_espacial_tmax_MOD11.png4.1.1.2 Mapa de la zona de estudio.
if(!require(tmap)) suppressMessages( suppressWarnings( install.packages("tmap") ))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))
chile_zona_a <- sf::st_read(file.path(ruta_datos,"shp/zona_a.shp"),quiet=TRUE)
chile_zona_a <- sf::st_as_sf(chile_zona_a)
tm_shape(chile_zona_a)+tm_polygons("NAME_1", legend.show = F) + tm_text("NAME_1")
4.2 Modelo de Regresión Lineal.
4.2.1 Preparar los datos.
Leer datos de TSS desde raster en grabado disco para la zona a.
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
modis_tss_stack_zona_a <- brick(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack_zona_a.tif"))Leer datos de Tmax en achivo RDS y convertir este data.frame en objeto sf para su uso posterior.
tmax_mod11_2017 <- readRDS(file.path(ruta_datos,"rds/data_unida_tmax_mod11_2017.rds"))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))
tmax_mod11_2017_sf <- st_as_sf( tmax_mod11_2017, coords=c('longitud','latitud'))
# tmax_mod11_2017_sf <- st_transform( tmax_mod11_2017_sf, CRS("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext"))4.2.2 Modelo de regresión lineal.
Análisis de los histogramas para ver si es necesario transformar las variables.
4.2.2.1 Modelos lineales para Enero.
4.2.2.1.1 Visualizar variable Tmax para Enero.
tmaxJan <- tmax_mod11_2017$Jan2017
id <- which(!is.na(tmaxJan))
tmaxJan <- tmaxJan[id]Histograma y densidad.
hist(tmaxJan)
plot(density(tmaxJan))
Comprobar normalidad.
qqnorm(tmaxJan)
shapiro.test(tmaxJan)##
## Shapiro-Wilk normality test
##
## data: tmaxJan
## W = 0.87335, p-value = 0.007432Creamos una segunda variable para Enero mediante transformación. Convertimos la variable Tmax utilizando sqrt(max(x+1) - x).
tmaxJant <- sqrt(max(tmaxJan+1)-tmaxJan)
hist(tmaxJant)
plot(density(tmaxJant))
qqnorm(tmaxJant)
shapiro.test(tmaxJant)##
## Shapiro-Wilk normality test
##
## data: tmaxJant
## W = 0.93613, p-value = 0.14844.2.2.1.2 Ajustar modelo de regresión para Enero.
tmaxJanMod <- tmax_mod11_2017$TmodisJan
tmaxJanMod <- tmaxJanMod[id]dataJan <- data.frame(tmaxEsta = tmaxJant,tmaxMod = tmaxJanMod)Modelo con valores transformados.
modJan <- lm(tmaxEsta~tmaxMod,dataJan)
summary(modJan)##
## Call:
## lm(formula = tmaxEsta ~ tmaxMod, data = dataJan)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4918 -0.4485 -0.2427 0.3153 1.3720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.43704 1.61598 3.365 0.00293 **
## tmaxMod -0.08619 0.05025 -1.715 0.10101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7572 on 21 degrees of freedom
## Multiple R-squared: 0.1229, Adjusted R-squared: 0.08112
## F-statistic: 2.942 on 1 and 21 DF, p-value: 0.101dataJan$fitted <- modJan$fitted.values
dataJan$res <- resid(modJan)
plot(dataJan$tmaxMod,dataJan$tmaxEsta)
Modelo con valores originales.
dataJan2 <- data.frame(tmaxEsta = tmaxJan,tmaxMod = tmaxJanMod)
modJan2 <- lm(tmaxEsta~tmaxJanMod,dataJan2)
summary(modJan2)##
## Call:
## lm(formula = tmaxEsta ~ tmaxJanMod, data = dataJan2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.806 -1.460 1.574 2.822 5.807
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.0985 9.3192 1.513 0.145
## tmaxJanMod 0.4453 0.2898 1.537 0.139
##
## Residual standard error: 4.367 on 21 degrees of freedom
## Multiple R-squared: 0.1011, Adjusted R-squared: 0.05828
## F-statistic: 2.362 on 1 and 21 DF, p-value: 0.1393dataJan2$fitted <- modJan2$fitted.values
dataJan2$res <- resid(modJan2)
plot(dataJan2$tmaxMod,dataJan2$tmaxEsta)
4.2.2.1.3 Dispersión de residuos para Enero.
if(!require(ggplot2)) suppressMessages( suppressWarnings( install.packages("ggplot2") ))ggplot(dataJan,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jan_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataJan,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jan_tss_vs_residuos.png"))## Saving 7 x 5 in imageggplot(dataJan2,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jan2_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataJan2,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jan2_tss_vs_residuos.png"))## Saving 7 x 5 in image4.2.2.2 Modelos lineales para Abril.
4.2.2.2.1 Visualizar variable Tmax para Abril.
tmaxApr <- tmax_mod11_2017$Apr2017
id <- which(!is.na(tmaxApr))
tmaxApr <- tmaxApr[id]Histograma y densidad.
hist(tmaxApr)
plot(density(tmaxApr))
Comprobar normalidad.
qqnorm(tmaxApr)
shapiro.test(tmaxApr)##
## Shapiro-Wilk normality test
##
## data: tmaxApr
## W = 0.80088, p-value = 0.0003958Creamos una segunda variable para Abril mediante transformación. Convertimos la variable Tmax utilizando sqrt(max(x+1) - x).
tmaxAprt <- sqrt(max(tmaxApr+1)-tmaxApr)
hist(tmaxAprt)
plot(density(tmaxAprt))
qqnorm(tmaxAprt)
shapiro.test(tmaxAprt)##
## Shapiro-Wilk normality test
##
## data: tmaxAprt
## W = 0.90298, p-value = 0.029134.2.2.2.2 Ajustar modelo de regresión para Abril.
tmaxAprMod <- tmax_mod11_2017$TmodisApr
tmaxAprMod <- tmaxAprMod[id]dataApr <- data.frame(tmaxEsta = tmaxAprt,tmaxMod = tmaxAprMod)Modelo con valores transformados.
modApr <- lm(tmaxEsta~tmaxMod,dataApr)
summary(modApr)##
## Call:
## lm(formula = tmaxEsta ~ tmaxMod, data = dataApr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82061 -0.26825 -0.06364 0.27448 0.78159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.83614 0.53384 10.932 3.98e-10 ***
## tmaxMod -0.20709 0.02956 -7.007 6.44e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4405 on 21 degrees of freedom
## Multiple R-squared: 0.7004, Adjusted R-squared: 0.6861
## F-statistic: 49.09 on 1 and 21 DF, p-value: 6.44e-07dataApr$fitted <- modApr$fitted.values
dataApr$res <- resid(modApr)
plot(dataApr$tmaxMod,dataApr$tmaxEsta)
Modelo con valores originales.
dataApr2 <- data.frame(tmaxEsta = tmaxApr,tmaxMod = tmaxAprMod)
modApr2 <- lm(tmaxEsta~tmaxAprMod,dataApr2)
summary(modApr2)##
## Call:
## lm(formula = tmaxEsta ~ tmaxAprMod, data = dataApr2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7486 -1.2989 0.1277 1.0643 3.6728
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6939 2.6512 0.262 0.796
## tmaxAprMod 1.0466 0.1468 7.130 4.95e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.188 on 21 degrees of freedom
## Multiple R-squared: 0.7077, Adjusted R-squared: 0.6937
## F-statistic: 50.84 on 1 and 21 DF, p-value: 4.954e-07dataApr2$fitted <- modApr2$fitted.values
dataApr2$res <- resid(modApr2)
plot(dataApr2$tmaxMod,dataApr2$tmaxEsta)
4.2.2.2.3 Dispersión de residuos para Abril.
if(!require(ggplot2)) suppressMessages( suppressWarnings( install.packages("ggplot2") ))ggplot(dataApr,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Apr_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataApr,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Apr_tss_vs_residuos.png"))## Saving 7 x 5 in imageggplot(dataApr2,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Apr2_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataApr2,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Apr2_tss_vs_residuos.png"))## Saving 7 x 5 in image4.2.2.3 Modelos lineales para Julio.
4.2.2.3.1 Visualizar variable Tmax para Julio.
tmaxJul <- tmax_mod11_2017$Jul2017
id <- which(!is.na(tmaxJul))
tmaxJul <- tmaxJul[id]Histograma y densidad.
hist(tmaxJul)
plot(density(tmaxJul))
Comprobar normalidad.
qqnorm(tmaxJul)
shapiro.test(tmaxJul)##
## Shapiro-Wilk normality test
##
## data: tmaxJul
## W = 0.86651, p-value = 0.005502Creamos una segunda variable para Julio mediante transformación. Convertimos la variable Tmax utilizando sqrt(max(x+1) - x).
tmaxJult <- sqrt(max(tmaxJul+1)-tmaxJul)
hist(tmaxJult)
plot(density(tmaxJult))
qqnorm(tmaxJult)
shapiro.test(tmaxJult)##
## Shapiro-Wilk normality test
##
## data: tmaxJult
## W = 0.95094, p-value = 0.30614.2.2.3.2 Ajustar modelo de regresión para Julio.
tmaxJulMod <- tmax_mod11_2017$TmodisJul
tmaxJulMod <- tmaxJulMod[id]dataJul <- data.frame(tmaxEsta = tmaxJult,tmaxMod = tmaxJulMod)Modelo con valores transformados.
modJul <- lm(tmaxEsta~tmaxMod,dataJul)
summary(modJul)##
## Call:
## lm(formula = tmaxEsta ~ tmaxMod, data = dataJul)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72594 -0.08366 0.03200 0.15567 0.63869
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.86096 0.11003 26.003 < 2e-16 ***
## tmaxMod -0.09827 0.01308 -7.515 2.21e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.342 on 21 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.716
## F-statistic: 56.47 on 1 and 21 DF, p-value: 2.211e-07dataJul$fitted <- modJul$fitted.values
dataJul$res <- resid(modJul)
plot(dataJul$tmaxMod,dataJul$tmaxEsta)
Modelo con valores originales.
dataJul2 <- data.frame(tmaxEsta = tmaxJul,tmaxMod = tmaxJulMod)
modJul2 <- lm(tmaxEsta~tmaxJulMod,dataJul2)
summary(modJul2)##
## Call:
## lm(formula = tmaxEsta ~ tmaxJulMod, data = dataJul2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7666 -0.6010 0.0251 0.5827 3.4581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.49584 0.49404 17.20 7.51e-14 ***
## tmaxJulMod 0.50262 0.05872 8.56 2.74e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.536 on 21 degrees of freedom
## Multiple R-squared: 0.7772, Adjusted R-squared: 0.7666
## F-statistic: 73.27 on 1 and 21 DF, p-value: 2.738e-08dataJul2$fitted <- modJul2$fitted.values
dataJul2$res <- resid(modJul2)
plot(dataJul2$tmaxMod,dataJul2$tmaxEsta)
4.2.2.3.3 Dispersión de residuos para Julio.
if(!require(ggplot2)) suppressMessages( suppressWarnings( install.packages("ggplot2") ))ggplot(dataJul,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jul_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataJul,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jul_tss_vs_residuos.png"))## Saving 7 x 5 in imageggplot(dataJul2,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jul2_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataJul2,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Jul2_tss_vs_residuos.png"))## Saving 7 x 5 in image4.2.2.4 Modelos lineales para Octubre.
4.2.2.4.1 Visualizar variable Tmax para Octubre.
tmaxOct <- tmax_mod11_2017$Oct2017
id <- which(!is.na(tmaxOct))
tmaxOct <- tmaxOct[id]Histograma y densidad.
hist(tmaxOct)
plot(density(tmaxOct))
Comprobar normalidad.
qqnorm(tmaxOct)
shapiro.test(tmaxOct)##
## Shapiro-Wilk normality test
##
## data: tmaxOct
## W = 0.78313, p-value = 0.0002083Creamos una segunda variable para Octubre mediante transformación. Convertimos la variable Tmax utilizando sqrt(max(x+1) - x).
tmaxOctt <- sqrt(max(tmaxOct+1)-tmaxOct)
hist(tmaxOctt)
plot(density(tmaxOctt))
qqnorm(tmaxOctt)
shapiro.test(tmaxOctt)##
## Shapiro-Wilk normality test
##
## data: tmaxOctt
## W = 0.88319, p-value = 0.011574.2.2.4.2 Ajustar modelo de regresión para Octubre.
tmaxOctMod <- tmax_mod11_2017$TmodisOct
tmaxOctMod <- tmaxOctMod[id]dataOct <- data.frame(tmaxEsta = tmaxOctt,tmaxMod = tmaxOctMod)Modelo con valores transformados.
modOct <- lm(tmaxEsta~tmaxMod,dataOct)
summary(modOct)##
## Call:
## lm(formula = tmaxEsta ~ tmaxMod, data = dataOct)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7704 -0.1644 0.0088 0.1595 0.6094
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.53257 0.19333 23.44 < 2e-16 ***
## tmaxMod -0.11891 0.01014 -11.73 1.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2857 on 21 degrees of freedom
## Multiple R-squared: 0.8675, Adjusted R-squared: 0.8612
## F-statistic: 137.5 on 1 and 21 DF, p-value: 1.113e-10dataOct$fitted <- modOct$fitted.values
dataOct$res <- resid(modOct)
plot(dataOct$tmaxMod,dataOct$tmaxEsta)
Modelo con valores originales.
dataOct2 <- data.frame(tmaxEsta = tmaxOct,tmaxMod = tmaxOctMod)
modOct2 <- lm(tmaxEsta~tmaxOctMod,dataOct2)
summary(modOct2)##
## Call:
## lm(formula = tmaxEsta ~ tmaxOctMod, data = dataOct2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7320 -0.6663 -0.2441 1.0254 2.2682
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.66211 0.92293 5.051 5.31e-05 ***
## tmaxOctMod 0.67172 0.04841 13.877 4.78e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.364 on 21 degrees of freedom
## Multiple R-squared: 0.9017, Adjusted R-squared: 0.897
## F-statistic: 192.6 on 1 and 21 DF, p-value: 4.78e-12dataOct2$fitted <- modOct2$fitted.values
dataOct2$res <- resid(modOct2)
plot(dataOct2$tmaxMod,dataOct2$tmaxEsta)
4.2.2.4.3 Dispersión de residuos para Octubre.
if(!require(ggplot2)) suppressMessages( suppressWarnings( install.packages("ggplot2") ))ggplot(dataOct,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Oct_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataOct,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Oct_tss_vs_residuos.png"))## Saving 7 x 5 in imageggplot(dataOct2,aes(tmaxMod,tmaxEsta)) +
geom_point() +
geom_smooth(method ='lm',formula=y~x) +
labs(x='TSS [°C]',y='Tmax [°C]') +
annotate("text",x=20,y=10,label="R^2==0.65",parse =TRUE) +
annotate("text",x=20,y=9.2,label="p-value < 0.01") +
annotate("text",x=12,y=25,label="Tmax = 3.3 + 0.92 TSS") +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Oct2_tss_vs_tmax_adjusted.png"))## Saving 7 x 5 in imageggplot(dataOct2,aes(tmaxMod,res)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(x='TSS [°C]',y='residuos') +
theme_bw()
ggsave(file.path(ruta_datos,"figs/Oct2_tss_vs_residuos.png"))## Saving 7 x 5 in image4.2.3 Predicción espacial no geoestadística.
4.2.3.1 Predicciones para Enero.
4.2.3.1.1 Variable predictora espacial para Enero.
Predicción espacial utilizando el modelo de regresión calculado.
modis_tss_stack_zona_a <- brick(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack_zona_a.tif"))
modis_tss_stack_zona_a_Jan <- subset(modis_tss_stack_zona_a,1)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
modis_tss_stack_zona_a_Jan <- modis_tss_stack_zona_a_Jan - 273.15plot(modis_tss_stack_zona_a_Jan)
names(modis_tss_stack_zona_a_Jan) <- 'tmaxMod'
modis_tss_stack_zona_a_Jan## class : RasterLayer
## dimensions : 54, 55, 2970 (nrow, ncol, ncell)
## resolution : 0.05, 0.05 (x, y)
## extent : -72.75, -70, -36.5, -33.8 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : tmaxMod
## values : 10.05001, 44.80999 (min, max)Vamos a utilizar predict para predecir en la grilla de modis.
Creamos máscara con valores de 1 dentro de los límites de las regiones y valores NA fuera de los límites.
mas <- !is.na(modis_tss_stack_zona_a_Jan)
mas[mas==0] <- NAPara poder aplicar el modelo de regresión en todo el dominio convertimos los NA en 0.
modis_tss_stack_zona_a_Jan[is.na(modis_tss_stack_zona_a_Jan)] <- 0
plot(modis_tss_stack_zona_a_Jan)
Convertimos el dominio que es objeto raster a objeto SpatialGriddataFrame ya que es el objeto que nos permite utilizar el paquete {gstat}.
grid_Jan <- as(modis_tss_stack_zona_a_Jan,'SpatialGridDataFrame')4.2.3.1.2 Predicción espacial para Enero.
Predicción espacial utilizando el modelo de regresión lineal almacenado en la variable modJan.
grid_Jan$pred <- predict(modJan,grid_Jan)
grid_Jan$se.fit <- predict(modJan,grid_Jan,se.fit = TRUE)$se.fitplot(raster(grid_Jan['pred']))
plot(raster(grid_Jan['se.fit']))
Raster stack con los valores de Tmax modis en la banda 1, valores de tmax estimada con el modelo de regresion para el mes de abril en la banda 2 y el error estandar en la banda 3.
tmaxEstimada <- stack(grid_Jan)
tmaxEstimada <- mask(tmaxEstimada,mas)Raster stack con tres bandas, la primera banda tiene la TSS de MODIS la segunda banda la Tmax estimada utilizando el modelo de regresion y la banda 3 tiene el error estandar.
plot(tmaxEstimada)
4.2.3.1.3 Predicción espacial utilizando IDW para Enero.
dataIDW <- tmax_mod11_2017_sf[,4] #datos espaciales para aplicar IDW
id <- which(!is.na(dataIDW$Jan2017)) #posición de los que no son NA
dataIDW <- dataIDW[id,] # dejo solo los valores que no son NAif(!require(gstat)) suppressMessages( suppressWarnings( install.packages("gstat") ))
grid_Jan$predIDW <- krige(Jan2017~1,as(dataIDW,'Spatial'),grid_Jan) #aplico IDW## [inverse distance weighted interpolation]Lo guardo en un objeto raster.
res_Jan <- stack(grid_Jan$predIDW)
res_Jan <- mask(res_Jan,mas)
plot(res_Jan[[1]])
Ploteo la diferencia entre los predicho de Tmax con el modelo de reggresión lineal y lo predicho con IDW
plot(tmaxEstimada$pred-res_Jan$var1.pred)
4.2.3.2 Predicciones para Abril.
4.2.3.2.1 Variable predictora espacial para Abril.
Predicción espacial utilizando el modelo de regresión calculado.
modis_tss_stack_zona_a <- brick(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack_zona_a.tif"))
modis_tss_stack_zona_a_Apr <- subset(modis_tss_stack_zona_a,2)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
modis_tss_stack_zona_a_Apr <- modis_tss_stack_zona_a_Apr - 273.15plot(modis_tss_stack_zona_a_Apr)
names(modis_tss_stack_zona_a_Apr) <- 'tmaxMod'
modis_tss_stack_zona_a_Apr## class : RasterLayer
## dimensions : 54, 55, 2970 (nrow, ncol, ncell)
## resolution : 0.05, 0.05 (x, y)
## extent : -72.75, -70, -36.5, -33.8 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : tmaxMod
## values : 3.309991, 28.55001 (min, max)Vamos a utilizar predict para predecir en la grilla de modis.
Creamos máscara con valores de 1 dentro de los límites de las regiones y valores NA fuera de los límites.
mas <- !is.na(modis_tss_stack_zona_a_Apr)
mas[mas==0] <- NAPara poder aplicar el modelo de regresión en todo el dominio convertimos los NA en 0.
modis_tss_stack_zona_a_Apr[is.na(modis_tss_stack_zona_a_Apr)] <- 0
plot(modis_tss_stack_zona_a_Apr)
Convertimos el dominio que es objeto raster a objeto SpatialGriddataFrame ya que es el objeto que nos permite utilizar el paquete {gstat}.
grid_Apr <- as(modis_tss_stack_zona_a_Apr,'SpatialGridDataFrame')4.2.3.2.2 Predicción espacial para Abril.
Predicción espacial utilizando el modelo de regresión lineal almacenado en la variable modApr.
grid_Apr$pred <- predict(modApr,grid_Apr)
grid_Apr$se.fit <- predict(modApr,grid_Apr,se.fit = TRUE)$se.fitplot(raster(grid_Apr['pred']))
plot(raster(grid_Apr['se.fit']))
Raster stack con los valores de Tmax modis en la banda 1, valores de tmax estimada con el modelo de regresion para el mes de abril en la banda 2 y el error estandar en la banda 3.
tmaxEstimada <- stack(grid_Apr)
tmaxEstimada <- mask(tmaxEstimada,mas)Raster stack con tres bandas, la primera banda tiene la TSS de MODIS la segunda banda la Tmax estimada utilizando el modelo de regresion y la banda 3 tiene el error estandar.
plot(tmaxEstimada)
4.2.3.2.3 Predicción espacial utilizando IDW para Abril.
dataIDW <- tmax_mod11_2017_sf[,5] #datos espaciales para aplicar IDW
id <- which(!is.na(dataIDW$Apr2017)) #posición de los que no son NA
dataIDW <- dataIDW[id,] # dejo solo los valores que no son NAif(!require(gstat)) suppressMessages( suppressWarnings( install.packages("gstat") ))
grid_Apr$predIDW <- krige(Apr2017~1,as(dataIDW,'Spatial'),grid_Apr) #aplico IDW## [inverse distance weighted interpolation]Lo guardo en un objeto raster.
res_Apr <- stack(grid_Apr$predIDW)
res_Apr <- mask(res_Apr,mas)
plot(res_Apr[[1]])
Ploteo la diferencia entre los predicho de Tmax con el modelo de reggresión lineal y lo predicho con IDW
plot(tmaxEstimada$pred-res_Apr$var1.pred)
4.2.3.3 Predicciones para Julio.
4.2.3.3.1 Variable predictora espacial para Julio.
Predicción espacial utilizando el modelo de regresión calculado.
modis_tss_stack_zona_a <- brick(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack_zona_a.tif"))
modis_tss_stack_zona_a_Jul <- subset(modis_tss_stack_zona_a,3)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
modis_tss_stack_zona_a_Jul <- modis_tss_stack_zona_a_Jul - 273.15plot(modis_tss_stack_zona_a_Jul)
names(modis_tss_stack_zona_a_Jul) <- 'tmaxMod'
modis_tss_stack_zona_a_Jul## class : RasterLayer
## dimensions : 54, 55, 2970 (nrow, ncol, ncell)
## resolution : 0.05, 0.05 (x, y)
## extent : -72.75, -70, -36.5, -33.8 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : tmaxMod
## values : -11.25001, 14.07 (min, max)Vamos a utilizar predict para predecir en la grilla de modis.
Creamos máscara con valores de 1 dentro de los límites de las regiones y valores NA fuera de los límites.
mas <- !is.na(modis_tss_stack_zona_a_Jul)
mas[mas==0] <- NAPara poder aplicar el modelo de regresión en todo el dominio convertimos los NA en 0.
modis_tss_stack_zona_a_Jul[is.na(modis_tss_stack_zona_a_Jul)] <- 0
plot(modis_tss_stack_zona_a_Jul)
Convertimos el dominio que es objeto raster a objeto SpatialGriddataFrame ya que es el objeto que nos permite utilizar el paquete {gstat}.
grid_Jul <- as(modis_tss_stack_zona_a_Jul,'SpatialGridDataFrame')4.2.3.3.2 Predicción espacial para Julio.
Predicción espacial utilizando el modelo de regresión lineal almacenado en la variable modJul.
grid_Jul$pred <- predict(modJul,grid_Jul)
grid_Jul$se.fit <- predict(modJul,grid_Jul,se.fit = TRUE)$se.fitplot(raster(grid_Jul['pred']))
plot(raster(grid_Jul['se.fit']))
Raster stack con los valores de Tmax modis en la banda 1, valores de tmax estimada con el modelo de regresion para el mes de abril en la banda 2 y el error estandar en la banda 3.
tmaxEstimada <- stack(grid_Jul)
tmaxEstimada <- mask(tmaxEstimada,mas)Raster stack con tres bandas, la primera banda tiene la TSS de MODIS la segunda banda la Tmax estimada utilizando el modelo de regresion y la banda 3 tiene el error estandar.
plot(tmaxEstimada)
4.2.3.3.3 Predicción espacial utilizando IDW para Julio.
dataIDW <- tmax_mod11_2017_sf[,6] #datos espaciales para aplicar IDW
id <- which(!is.na(dataIDW$Jul2017)) #posición de los que no son NA
dataIDW <- dataIDW[id,] # dejo solo los valores que no son NAif(!require(gstat)) suppressMessages( suppressWarnings( install.packages("gstat") ))
grid_Jul$predIDW <- krige(Jul2017~1,as(dataIDW,'Spatial'),grid_Jul) #aplico IDW## [inverse distance weighted interpolation]Lo guardo en un objeto raster.
res_Jul <- stack(grid_Jul$predIDW)
res_Jul <- mask(res_Jul,mas)
plot(res_Jul[[1]])
Ploteo la diferencia entre los predicho de Tmax con el modelo de reggresión lineal y lo predicho con IDW
plot(tmaxEstimada$pred-res_Jul$var1.pred)
4.2.3.4 Predicciones para Octubre.
4.2.3.4.1 Variable predictora espacial para Octubre.
Predicción espacial utilizando el modelo de regresión calculado.
modis_tss_stack_zona_a <- brick(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack_zona_a.tif"))
modis_tss_stack_zona_a_Oct <- subset(modis_tss_stack_zona_a,4)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
modis_tss_stack_zona_a_Oct <- modis_tss_stack_zona_a_Oct - 273.15plot(modis_tss_stack_zona_a_Oct)
names(modis_tss_stack_zona_a_Oct) <- 'tmaxMod'
modis_tss_stack_zona_a_Oct## class : RasterLayer
## dimensions : 54, 55, 2970 (nrow, ncol, ncell)
## resolution : 0.05, 0.05 (x, y)
## extent : -72.75, -70, -36.5, -33.8 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : tmaxMod
## values : -1.769995, 28.70999 (min, max)Vamos a utilizar predict para predecir en la grilla de modis.
Creamos máscara con valores de 1 dentro de los límites de las regiones y valores NA fuera de los límites.
mas <- !is.na(modis_tss_stack_zona_a_Oct)
mas[mas==0] <- NAPara poder aplicar el modelo de regresión en todo el dominio convertimos los NA en 0.
modis_tss_stack_zona_a_Oct[is.na(modis_tss_stack_zona_a_Oct)] <- 0
plot(modis_tss_stack_zona_a_Oct)
Convertimos el dominio que es objeto raster a objeto SpatialGriddataFrame ya que es el objeto que nos permite utilizar el paquete {gstat}.
grid_Oct <- as(modis_tss_stack_zona_a_Oct,'SpatialGridDataFrame')4.2.3.4.2 Predicción espacial para Octubre.
Predicción espacial utilizando el modelo de regresión lineal almacenado en la variable modOct.
grid_Oct$pred <- predict(modOct,grid_Oct)
grid_Oct$se.fit <- predict(modOct,grid_Oct,se.fit = TRUE)$se.fitplot(raster(grid_Oct['pred']))
plot(raster(grid_Oct['se.fit']))
Raster stack con los valores de Tmax modis en la banda 1, valores de tmax estimada con el modelo de regresion para el mes de abril en la banda 2 y el error estandar en la banda 3.
tmaxEstimada <- stack(grid_Oct)
tmaxEstimada <- mask(tmaxEstimada,mas)Raster stack con tres bandas, la primera banda tiene la TSS de MODIS la segunda banda la Tmax estimada utilizando el modelo de regresion y la banda 3 tiene el error estandar.
plot(tmaxEstimada)
4.2.3.4.3 Predicción espacial utilizando IDW para Octubre.
dataIDW <- tmax_mod11_2017_sf[,7] #datos espaciales para aplicar IDW
id <- which(!is.na(dataIDW$Oct2017)) #posición de los que no son NA
dataIDW <- dataIDW[id,] # dejo solo los valores que no son NAif(!require(gstat)) suppressMessages( suppressWarnings( install.packages("gstat") ))
grid_Oct$predIDW <- krige(Oct2017~1,as(dataIDW,'Spatial'),grid_Oct) #aplico IDW## [inverse distance weighted interpolation]Lo guardo en un objeto raster.
res_Oct <- stack(grid_Oct$predIDW)
res_Oct <- mask(res_Oct,mas)
plot(res_Oct[[1]])
Ploteo la diferencia entre los predicho de Tmax con el modelo de reggresión lineal y lo predicho con IDW
plot(tmaxEstimada$pred-res_Oct$var1.pred)
4.2.4 Resultados del Taller 2.
4.2.4.1 Gráficas de variación de residuos versus TSS.
png(file.path(ruta_datos,"figs/grafico_residuos_tss.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(dataJan2$tmaxMod,dataJan2$res,main="Residuos vs TSS Jan",xlab="TSS",ylab="Residuos")
abline(h=0,col="red")
plot(dataApr2$tmaxMod,dataApr2$res,main="Residuos vs TSS Apr",xlab="TSS",ylab="Residuos")
abline(h=0,col="red")
plot(dataJul2$tmaxMod,dataJul2$res,main="Residuos vs TSS Jul",xlab="TSS",ylab="Residuos")
abline(h=0,col="red")
plot(dataOct2$tmaxMod,dataOct2$res,main="Residuos vs TSS Oct",xlab="TSS",ylab="Residuos")
abline(h=0,col="red")
suppressMessages( dev.off())## png
## 2
4.2.4.2 Gráficas de ajuste de la regresión.
if(!require(visreg)) suppressMessages( suppressWarnings( install.packages("visreg") ))## Loading required package: visregif(!require(broom)) suppressMessages( suppressWarnings( install.packages("broom") ))## Loading required package: broompng(file.path(ruta_datos,"figs/grafico_ols_intervalos.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
visreg(modJan2,main="OLS Jan",xlab="TSS",ylab="Tmax")
legend('topleft',legend=c(paste("R2:",round(summary(modJan2)$r.squared,4)),
paste("Pval:",round(glance(modJan2)$p.value,4))),text.col=c("red","blue"),cex=.8)
visreg(modApr2,main="OLS Apr",xlab="TSS",ylab="Tmax")
legend('topleft',legend=c(paste("R2:",round(summary(modApr2)$r.squared,4)),
paste("Pval:",round(glance(modApr2)$p.value,4))),text.col=c("red","blue"),cex=.8)
visreg(modJul2,main="OLS Jul",xlab="TSS",ylab="Tmax")
legend('topleft',legend=c(paste("R2:",round(summary(modJul2)$r.squared,4)),
paste("Pval:",round(glance(modJul2)$p.value,4))),text.col=c("red","blue"),cex=.8)
visreg(modOct2,main="OLS Oct",xlab="TSS",ylab="Tmax")
legend('topleft',legend=c(paste("R2:",round(summary(modOct2)$r.squared,4)),
paste("Pval:",round(glance(modOct2)$p.value,4))),text.col=c("red","blue"),cex=.8)
suppressMessages( dev.off())## png
## 2
4.2.4.3 Mapas de predicción espacial.
png(file.path(ruta_datos,"figs/mapa_prediccion_ols.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(raster(grid_Jan['pred']),main="Pred OLS Jan")
plot(raster(grid_Apr['pred']),main="Pred OLS Apr")
plot(raster(grid_Jul['pred']),main="Pred OLS Jul")
plot(raster(grid_Oct['pred']),main="Pred OLS Oct")
suppressMessages( dev.off())## png
## 2
png(file.path(ruta_datos,"figs/mapa_prediccion_idw.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(res_Jan[[1]],main="Pred IDW Jan")
plot(res_Apr[[1]],main="Pred IDW Apr")
plot(res_Jul[[1]],main="Pred IDW Jul")
plot(res_Oct[[1]],main="Pred IDW Oct")
suppressMessages( dev.off())## png
## 2
4.3 Análisis de Variograma.
4.3.1 Preparar los datos.
Crear los data.frame con los datos de temperatura maxima observada y los valores extraidos de las variables predictoras TSS, distancia a la costa y elevación.
data_unida <- readRDS(file.path(ruta_datos,"rds/data_unida_tmax_mod11_2017.rds"))if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))4.3.1.1 Raster de distancia a la costa.
Creamos raster de distancia a la costa que se usará como predictor de temperatura máxima mensual.
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))Descargamos polígono de Chile, level=0 indica que la división es a nivel país, sin regiones ni comunas.
if(file.exists(file.path(ruta_datos,"rds/gadm36_CHL_0_sp.rds")))
{ chl <- readRDS(file.path(ruta_datos,"rds/gadm36_CHL_0_sp.rds"))
} else {
chl <- getData('GADM',country='chl',level=0)
}Convertimos a objeto sf.
chl_sf <- st_as_sf(chl)
class(chl_sf)## [1] "sf" "data.frame"Leemos un archivo box_costa que contiene los límites de la zona a la cual queremos extraer la línea de costa.
box_costa_sf <- st_read(file.path(ruta_datos,"shp/box_costa.shp"),quiet=TRUE)
class(box_costa_sf)## [1] "sf" "data.frame"Graficamos las geometrías de los objetos sf.
plot(chl_sf$geometry)
plot(box_costa_sf$geometry,add=TRUE,col='red')
Intersección de borde de chile con box de área de costa convertimos de MULIPOLIGONO a MULTILINESTRING.
chl_sf_line <- st_cast(chl_sf,'MULTILINESTRING')Intersectamos las líneas del borde de Chile con el box.
borde_sf <- st_intersection(chl_sf_line,box_costa_sf)## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all geometriesDescargamos las regiones de Chile en nuestra zona de estudio, con la opción level=1, obtenemos los límites regionales.
if(file.exists(file.path(ruta_datos,"rds/gadm36_CHL_1_sp.rds")))
{ regs <- readRDS(file.path(ruta_datos,"rds/gadm36_CHL_1_sp.rds"))
} else {
regs <- getData('GADM',country='chl',level=1)
}Convertimos en sf.
regs_sf <- st_as_sf(regs)Seleccionamos BioBio.
regs_sf_sel <- regs_sf[c(6),c(1)]Seleccionamos la zona a.
zona_a <- readRDS(file.path(ruta_datos,"rds/chile_zona_a.rds"))
zona_a <- st_union(zona_a)## although coordinates are longitude/latitude, st_union assumes that they are planarIntersectamos borde de costa con la región seleccionada.
bor_reg <- st_intersection(borde_sf,zona_a)## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all geometries# bor_reg <- st_transform(bor_reg, CRS("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext"))Unimos todas las líneas en una sola de bor_reg.
bor_reg <- st_union(bor_reg)## although coordinates are longitude/latitude, st_union assumes that they are planarObjeto raster que vamos a utilizar para calcular las distancias y guardarlas.
modis_tss_stack <- raster(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack.tif"))Lo cortamos y aplicamos máscara para la región.
extent(modis_tss_stack)## class : Extent
## xmin : -180
## xmax : 180
## ymin : -90
## ymax : 90extent(as(zona_a,'Spatial'))## class : Extent
## xmin : -72.77348
## xmax : -70.01145
## ymin : -36.48368
## ymax : -33.80633crs(modis_tss_stack)## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0crs(as(zona_a,'Spatial'))## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0modis_tss_stack <- crop(modis_tss_stack,as(zona_a,'Spatial'))
modis_tss_stack <- mask(modis_tss_stack,as(zona_a,'Spatial'))Creamos un raster para el cálculo de distancia.
distRas <- modis_tss_stackCambiamos los valores para que el raster tenga valores únicos en cada pixel.
distRas[] <- 1:ncell(distRas)Conviertimos de raster a vector, objeto SpatialPolygon.
polDist <- rasterToPolygons(distRas)Lo conviertimos en objeto sf.
polDist <- st_as_sf(polDist)
plot(polDist$geometry)
Calculamos los centroides de cada polígono.
centro <- st_centroid(polDist)## Warning in st_centroid.sf(polDist): st_centroid assumes attributes are constant over geometries of x## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon): st_centroid does not give correct centroids for longitude/latitude dataplot(centro$geometry)
Calculamos la distancia entre el borde de la costa y el centroide de cada polígono que representa a cada pixel de MOD11. Optamos por ejecutar la función st_distance sólo de forma condicional, puesto que tarda mucho tiempo. Así que la ejecutamos solamente la primera vez, en caso de no existir el objeto rds que la almacena, en caso contrario que la cargue en sesión desde dicho archivo.
if(file.exists(file.path(ruta_datos,"rds/dist_bor_reg_centro.rds")))
{ dist_bor_reg_centro <- readRDS(file.path(ruta_datos,"rds/dist_bor_reg_centro.rds"))
} else { dist_bor_reg_centro <- st_distance(bor_reg,centro)
saveRDS( dist_bor_reg_centro , file.path(ruta_datos,"rds/dist_bor_reg_centro.rds"))
}
plot(dist_bor_reg_centro)
# crs(bor_reg)
# crs(centro)saveRDS(bor_reg,file.path(ruta_datos,"rds/bor_reg_zona_b.rds"))saveRDS(centro,file.path(ruta_datos,"rds/centro_zona_b.rds"))Asignamos la distancia calculada a los pixeles del raster distRas.
distRas[] <- as.numeric(dist_bor_reg_centro)plot(distRas)
plot(st_geometry(bor_reg),add=TRUE)
plot(st_geometry(zona_a),add=TRUE)
Como usaremos las regiones seleccionadas para intersectarlas con la imagen modis, usaremos la proyección 4326, ya que la imagen modis también la grabamos en archivo en disco con esa proyección.
# zona_a <- st_transform( zona_a, CRS("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext"))
zona_a1 <- as( zona_a, 'Spatial' )
# zona_a1 <- spTransform( zona_a1, CRS(proj4string( distRas)))
distCosta <- crop(distRas , zona_a1)
distCosta <- mask(distCosta, zona_a1)
plot(distCosta)
Guardamos el raster de la distancia a la costa para poder usarlo más adelante.
writeRaster(distCosta,file.path(ruta_datos,"raster/dist_costa_biobio.tif"), overwrite=TRUE)4.3.1.2 Raster de elevación.
Descargamos los polígonos de las regiones de Chile.
if(file.exists(file.path(ruta_datos,"rds/gadm36_CHL_1_sp.rds")))
{ regs <- readRDS(file.path(ruta_datos,"rds/gadm36_CHL_1_sp.rds"))
} else {
regs <- getData('GADM',country='chl',level=1)
}
regs <- st_as_sf(regs)Seleccionamos la región del BioBio.
# zona_a <- regs_sf[c(6),c(1)]
st_centroid(zona_a)## Warning in st_centroid.sfc(zona_a): st_centroid does not give correct centroids for longitude/latitude data## Geometry set for 1 feature
## geometry type: POINT
## dimension: XY
## bbox: xmin: -71.29425 ymin: -35.17712 xmax: -71.29425 ymax: -35.17712
## CRS: 4326## POINT (-71.29425 -35.17712)Descargamos las elevaciones.
if(file.exists(file.path(ruta_datos,"grd/CHL1_msk_alt.grd")))
{ ele <- raster::raster(file.path(ruta_datos,"grd/CHL1_msk_alt.grd"))
} else {
ele <- getData('alt',country = 'CHL')
}
ele <- ele[[1]]Cortamos las elevaciones para la región elegida.
eleR <- crop(ele,as(zona_a,'Spatial'))
eleR <- mask(eleR,as(zona_a,'Spatial'))Cargamos el raster mod11 de generado anteriormente.
modis_tss_stack <- stack(file.path(ruta_datos,"raster/MOD11C3_A2017_tss_stack.tif"))
# crs(modis_tss_stack)Lo cortamos para la región de interés.
modis_tss_stackR <- crop(modis_tss_stack,as(zona_a,'Spatial'))
modis_tss_stackR <- mask(modis_tss_stackR,as(zona_a,'Spatial'))res(modis_tss_stackR)## [1] 0.05 0.05res(eleR)## [1] 0.008333333 0.008333333Hacemos un resampleo de la elevación para que tenga el mismo origen y resolución que el raster mod11.
eleRMOD <- resample(eleR,modis_tss_stackR)
plot(eleRMOD)
plot(modis_tss_stackR)
Ahora guardamos el raste de elevación creado para el área de estudio para ser usado más tarde.
writeRaster(eleRMOD,file.path(ruta_datos,"raster/elevacion_biobio.tif"), overwrite=TRUE)4.3.1.3 Usamos los rasters de distancia y elevación.
dist_costa_biobio <- raster(file.path(ruta_datos,"raster/dist_costa_biobio.tif"))
elevacion_biobio <- raster(file.path(ruta_datos,"raster/elevacion_biobio.tif"))pred_dist_ele_biobio <- brick(dist_costa_biobio,elevacion_biobio)dataEx <- raster::extract(pred_dist_ele_biobio,data_unida)
plot(pred_dist_ele_biobio)
plot(data_unida,add=TRUE)## Warning in plot.sf(data_unida, add = TRUE): ignoring all but the first attribute
En el código anterior nos encontramos con que no existen estaciones meteorológicas en la región de BioBio. Por este motivo realizaremos un recuento de cuántas estaciones tenemos por región en todo chile para decidir qué zona escoger para el estudio.
if(file.exists(file.path(ruta_datos,"rds/estaciones_chile.rds")))
{ estaciones_chile <- readRDS(file.path(ruta_datos,"rds/estaciones_chile.rds"))
} else {
tmax_2017 <- readRDS(file.path(ruta_datos,"rds/tmax_2017_taller1.rds"))
tmax_2017_sf <- sf::st_as_sf( tmax_2017, coords=c('longitud','latitud'))
st_crs(tmax_2017_sf) <- 4326
estaciones_chile <- sf::st_intersection(tmax_2017_sf,regs)
saveRDS(estaciones_chile,file.path(ruta_datos,"rds/estaciones_chile.rds"))
}Para contabilizar datos faltantes y no faltantes por región usaremos el formato longer que nos permita filtrar esa información más fácilmente.
if(!require(tidyverse)) suppressMessages( suppressWarnings( install.packages("tidyverse") ))## Loading required package: tidyverse## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✔ tibble 3.0.3 ✔ stringr 1.4.0
## ✔ readr 1.4.0 ✔ forcats 0.5.0
## ✔ purrr 0.3.4## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()estaciones_chile_nocoord <- st_drop_geometry(estaciones_chile)
estaciones_chile_longer <- pivot_longer(estaciones_chile_nocoord, 4:7,names_to='Meses',values_to = 'Tmax')
str(estaciones_chile_longer)## tibble [1,556 × 15] (S3: tbl_df/tbl/data.frame)
## $ codigo_estacion: chr [1:1556] "11043001" "11043001" "11043001" "11043001" ...
## $ institucion : chr [1:1556] "DGA" "DGA" "DGA" "DGA" ...
## $ altura : chr [1:1556] "10" "10" "10" "10" ...
## $ GID_0 : chr [1:1556] "CHL" "CHL" "CHL" "CHL" ...
## $ NAME_0 : chr [1:1556] "Chile" "Chile" "Chile" "Chile" ...
## $ GID_1 : chr [1:1556] "CHL.1_1" "CHL.1_1" "CHL.1_1" "CHL.1_1" ...
## $ NAME_1 : chr [1:1556] "Aisén del General Carlos Ibáñez del Campo" "Aisén del General Carlos Ibáñez del Campo" "Aisén del General Carlos Ibáñez del Campo" "Aisén del General Carlos Ibáñez del Campo" ...
## $ VARNAME_1 : chr [1:1556] "Aisén|Aysén|Aysén del General Carlos Ibáñez del Campo" "Aisén|Aysén|Aysén del General Carlos Ibáñez del Campo" "Aisén|Aysén|Aysén del General Carlos Ibáñez del Campo" "Aisén|Aysén|Aysén del General Carlos Ibáñez del Campo" ...
## $ NL_NAME_1 : chr [1:1556] NA NA NA NA ...
## $ TYPE_1 : chr [1:1556] "Región" "Región" "Región" "Región" ...
## $ ENGTYPE_1 : chr [1:1556] "Region" "Region" "Region" "Region" ...
## $ CC_1 : chr [1:1556] "XI" "XI" "XI" "XI" ...
## $ HASC_1 : chr [1:1556] "CL.AI" "CL.AI" "CL.AI" "CL.AI" ...
## $ Meses : chr [1:1556] "Jan2017" "Apr2017" "Jul2017" "Oct2017" ...
## $ Tmax : num [1:1556] NA NA NA NA 18 ...Contabilizamos las estaciones por región. Esto no descuenta aquellas estaciones sin datos faltantes.
estaciones_region <- table(estaciones_chile_longer$NAME_1,
estaciones_chile_longer$Meses)
estaciones_region <- as.data.frame.matrix( estaciones_region)
estaciones_region <- estaciones_region[,c(2,1,3,4)]
estaciones_region## Jan2017 Apr2017 Jul2017 Oct2017
## Aisén del General Carlos Ibáñez del Campo 42 42 42 42
## Antofagasta 38 38 38 38
## Araucanía 29 29 29 29
## Arica y Parinacota 23 23 23 23
## Atacama 23 23 23 23
## Bío-Bío 11 11 11 11
## Coquimbo 37 37 37 37
## Libertador General Bernardo O'Higgins 12 12 12 12
## Los Lagos 20 20 20 20
## Los Ríos 7 7 7 7
## Magallanes y Antártica Chilena 45 45 45 45
## Maule 25 25 25 25
## Ñuble 9 9 9 9
## Región Metropolitana de Santiago 33 33 33 33
## Tarapacá 14 14 14 14
## Valparaíso 21 21 21 21# dim( estaciones_chile)
# estaciones_chile [estaciones_chile$NAME_1=="Bío-Bío",]
# estaciones_chile [estaciones_chile$NAME_1=="Maule",]
# estaciones_chile [estaciones_chile$NAME_1=="Libertador General Bernardo O'Higgins",]Ahora filtramos las observaciones de Tmax faltantes.
estaciones_chile_longer_filter <- estaciones_chile_longer %>%
filter(!is.na(Tmax))
dim(estaciones_chile_longer_filter)## [1] 738 15tmax_region <- table(estaciones_chile_longer_filter$NAME_1,
estaciones_chile_longer_filter$Meses)
tmax_region <- as.data.frame.matrix( tmax_region)
tmax_region <- tmax_region[,c(2,1,3,4)]
tmax_region## Jan2017 Apr2017 Jul2017 Oct2017
## Aisén del General Carlos Ibáñez del Campo 18 17 19 20
## Antofagasta 9 9 9 11
## Araucanía 8 9 7 7
## Arica y Parinacota 8 8 7 8
## Atacama 10 9 7 8
## Bío-Bío 8 8 8 8
## Coquimbo 10 9 10 10
## Libertador General Bernardo O'Higgins 7 7 7 7
## Los Lagos 13 13 12 13
## Los Ríos 3 2 3 3
## Magallanes y Antártica Chilena 37 37 37 34
## Maule 16 16 16 16
## Ñuble 3 4 4 5
## Región Metropolitana de Santiago 19 20 19 19
## Tarapacá 4 4 4 4
## Valparaíso 13 13 13 124.3.1.4 Intersectar los datos TSS y Tmax.
dataJan <- cbind(data_unida[,c(4,8)],dataEx)
names(dataJan)[1:4] <- c('Tmax','TSS','dist','ele')
saveRDS(dataJan,file.path(ruta_datos,'rds/data_january.rds'))dataApr <- cbind(data_unida[,c(5,9)],dataEx)
names(dataApr)[1:4] <- c('Tmax','TSS','dist','ele')
saveRDS(dataApr,file.path(ruta_datos,'rds/data_april.rds'))dataJul <- cbind(data_unida[,c(6,10)],dataEx)
names(dataJul)[1:4] <- c('Tmax','TSS','dist','ele')
saveRDS(dataJul,file.path(ruta_datos,'rds/data_july.rds'))dataOct <- cbind(data_unida[,c(7,11)],dataEx)
names(dataOct)[1:4] <- c('Tmax','TSS','dist','ele')
saveRDS(dataOct,file.path(ruta_datos,'rds/data_october.rds'))4.3.2 Modelos de variograma.
4.3.2.1 Modelos de regresión múltiple.
modJan <- lm(Tmax~TSS+dist+ele,dataJan)
summary(modJan)##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataJan)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9884 -0.5993 0.0650 0.9315 2.3007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.130e+01 3.761e+00 5.664 2.26e-05 ***
## TSS 2.495e-01 1.144e-01 2.181 0.04273 *
## dist 5.320e-05 1.427e-05 3.728 0.00154 **
## ele -6.103e-03 6.428e-04 -9.494 1.97e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.564 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8968, Adjusted R-squared: 0.8796
## F-statistic: 52.14 on 3 and 18 DF, p-value: 4.456e-09modApr <- lm(Tmax~TSS+dist+ele,dataApr)
summary(modApr)##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataApr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5353 -0.4794 0.0719 0.9757 2.5336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.546e+01 4.408e+00 3.507 0.00251 **
## TSS 3.553e-01 1.989e-01 1.786 0.09089 .
## dist 9.152e-06 1.498e-05 0.611 0.54891
## ele -3.206e-03 8.497e-04 -3.774 0.00139 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.651 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8568, Adjusted R-squared: 0.8329
## F-statistic: 35.89 on 3 and 18 DF, p-value: 8.35e-08modJul <- lm(Tmax~TSS+dist+ele,dataJul)
summary(modJul)##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataJul)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6094 -0.6500 -0.0589 0.5406 4.0711
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.068e+00 1.942e+00 4.153 0.000597 ***
## TSS 5.809e-01 1.359e-01 4.276 0.000455 ***
## dist -9.957e-06 1.457e-05 -0.683 0.503020
## ele 9.247e-04 9.855e-04 0.938 0.360528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.608 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.7781, Adjusted R-squared: 0.7411
## F-statistic: 21.04 on 3 and 18 DF, p-value: 4.111e-06modOct <- lm(Tmax~TSS+dist+ele,dataOct)
summary(modOct)##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataOct)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.98716 -0.85230 0.09615 0.79237 1.59648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.711e+00 1.819e+00 5.338 4.50e-05 ***
## TSS 4.928e-01 7.762e-02 6.349 5.56e-06 ***
## dist -3.001e-06 1.135e-05 -0.264 0.7944
## ele -1.323e-03 7.294e-04 -1.813 0.0865 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.189 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9359, Adjusted R-squared: 0.9252
## F-statistic: 87.62 on 3 and 18 DF, p-value: 6.239e-114.3.2.2 Índice de Moran.
Instalamos el paquete spdep para obtener las funciones, moran.test y moran.plot.
if(!require(spdep)) suppressMessages( suppressWarnings( install.packages("spdep") ))4.3.2.2.1 Índice de Moran para Enero.
#moran.test(dataJan$Tmax,nsim=999,4.3.2.2.2 Índice de Moran para Abril.
#moran.test(dataJan$Tmax,nsim=999,4.3.2.2.3 Índice de Moran para Julio.
#moran.test(dataJan$Tmax,nsim=999,4.3.2.2.4 Índice de Moran para Octubre.
#moran.test(dataJan$Tmax,nsim=999,4.3.2.3 Variograma experimental usando distancia a la costa.
if(!require(gstat)) suppressMessages( suppressWarnings( install.packages("gstat") ))
options(warn = -1) # don't print warningsdataJan <- readRDS(file.path(ruta_datos,"rds/data_january.rds"))
dataApr <- readRDS(file.path(ruta_datos,"rds/data_april.rds"))
dataJul <- readRDS(file.path(ruta_datos,"rds/data_july.rds"))
dataOct <- readRDS(file.path(ruta_datos,"rds/data_october.rds"))
dataJan <- dataJan [ !is.na(dataJan$dist), ]
dataApr <- dataApr [ !is.na(dataApr$dist), ]
dataJul <- dataJul [ !is.na(dataJul$dist), ]
dataOct <- dataOct [ !is.na(dataOct$dist), ]4.3.2.3.1 Variograma experimental usando Tmax para Enero.
v.Jan <- variogram(Tmax~1,dataJan)
plot(v.Jan,plot.numbers=TRUE)
4.3.2.3.2 Variograma experimental usando Tmax para Abril.
v.Apr <- variogram(Tmax~1,dataApr)
plot(v.Apr,plot.numbers=TRUE)
4.3.2.3.3 Variograma experimental usando Tmax para Julio.
v.Jul <- variogram(Tmax~1,dataJul)
plot(v.Jul,plot.numbers=TRUE)
4.3.2.3.4 Variograma experimental usando Tmax para Octubre.
v.Oct <- variogram(Tmax~1,dataOct)
plot(v.Oct,plot.numbers=TRUE)
4.3.2.4 Variograma experimental usando distancia.
4.3.2.4.1 Variograma experimental usando distancia para Enero.
v.dist.Jan <- variogram(dist~1,dataJan)
plot(v.dist.Jan,plot.numbers=TRUE)
4.3.2.4.2 Variograma experimental usando distancia para Abril.
v.dist.Apr <- variogram(dist~1,dataApr)
plot(v.dist.Apr,plot.numbers=TRUE)
4.3.2.4.3 Variograma experimental usando distancia para Julio.
v.dist.Jul <- variogram(dist~1,dataJul)
plot(v.dist.Jul,plot.numbers=TRUE)
4.3.2.4.4 Variograma experimental usando distancia para Octubre.
v.dist.Oct <- variogram(dist~1,dataOct)
plot(v.dist.Oct,plot.numbers=TRUE)
4.3.2.5 Variograma experimental usando elevación.
4.3.2.5.1 Variograma experimental usando elevación para Enero.
v.ele.Jan <- variogram(ele~1,dataJan)
plot(v.ele.Jan,plot.numbers=TRUE)
4.3.2.5.2 Variograma experimental usando elevación para Abril.
v.ele.Apr <- variogram(ele~1,dataApr)
plot(v.ele.Apr,plot.numbers=TRUE)
4.3.2.5.3 Variograma experimental usando elevación para Julio.
v.ele.Jul <- variogram(ele~1,dataJul)
plot(v.ele.Jul,plot.numbers=TRUE)
4.3.2.5.4 Variograma experimental usando elevación para Octubre.
v.ele.Oct <- variogram(ele~1,dataOct)
plot(v.ele.Oct,plot.numbers=TRUE)
4.3.2.6 Variograma ajustado usando distancia a la costa.
4.3.2.6.1 Variograma ajustado usando distancia para Enero.
v.dist.Jan <- variogram(dist~1,dataJan)
v.dist.Jan.fit <- suppressWarnings( fit.variogram(v.dist.Jan, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.dist.Jan.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.6.2 Variograma ajustado usando distancia para Abril.
v.dist.Apr <- variogram(dist~1,dataApr)
v.dist.Apr.fit <- suppressWarnings( fit.variogram(v.dist.Apr, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.dist.Apr.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.6.3 Variograma ajustado usando distancia para Julio.
v.dist.Jul <- variogram(dist~1,dataJul)
v.dist.Jul.fit <- suppressWarnings( fit.variogram(v.dist.Jul, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.dist.Jul.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.6.4 Variograma ajustado usando distancia para Octubre.
v.dist.Oct <- variogram(dist~1,dataOct)
v.dist.Oct.fit <- suppressWarnings( fit.variogram(v.dist.Oct, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.dist.Oct.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.7 Variograma ajustado usando elevación.
4.3.2.7.1 Variograma ajustado usando elevación para Enero.
v.ele.Jan <- variogram(ele~1,dataJan)
v.ele.Jan.fit <- suppressWarnings( fit.variogram(v.ele.Jan, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.ele.Jan.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.7.2 Variograma ajustado usando elevación para Abril.
v.ele.Apr <- variogram(ele~1,dataApr)
v.ele.Apr.fit <- suppressWarnings( fit.variogram(v.ele.Apr, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.ele.Apr.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.7.3 Variograma ajustado usando elevación para Julio.
v.ele.Jul <- variogram(ele~1,dataJul)
v.ele.Jul.fit <- suppressWarnings( fit.variogram(v.ele.Jul, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.ele.Jul.fit,plot.numbers=TRUE,cutoff=10000)
4.3.2.7.4 Variograma ajustado usando elevación para Octubre.
v.ele.Oct <- variogram(ele~1,dataOct)
v.ele.Oct.fit <- suppressWarnings( fit.variogram(v.ele.Oct, vgm(c("Ste")), fit.kappa = seq(.3,5,.01)))
plot(v.ele.Oct.fit,plot.numbers=TRUE,cutoff=10000)
4.3.3 Resultados del taller 3.
4.3.3.1 Gráficos del índice de Moran.
4.3.3.2 Gráficos de variograma experimental usando distancia y elevación.
png(file.path(ruta_datos,"figs/variograma_experimental_tmax.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(v.Jan,plot.numbers=TRUE,cutoff=5000)
plot(v.Apr,plot.numbers=TRUE,cutoff=5000)
plot(v.Jul,plot.numbers=TRUE,cutoff=5000)
plot(v.Oct,plot.numbers=TRUE,cutoff=5000)
dev.off()## png
## 2
png(file.path(ruta_datos,"figs/variograma_experimental_distancia.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(v.dist.Jan,plot.numbers=TRUE,cutoff=5000)
plot(v.dist.Apr,plot.numbers=TRUE,cutoff=5000)
plot(v.dist.Jul,plot.numbers=TRUE,cutoff=5000)
plot(v.dist.Oct,plot.numbers=TRUE,cutoff=5000)
dev.off()## png
## 2
png(file.path(ruta_datos,"figs/variograma_experimental_elevacion.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(v.ele.Jan,plot.numbers=TRUE,cutoff=5000)
plot(v.ele.Apr,plot.numbers=TRUE,cutoff=5000)
plot(v.ele.Jul,plot.numbers=TRUE,cutoff=5000)
plot(v.ele.Oct,plot.numbers=TRUE,cutoff=5000)
dev.off()## png
## 2
4.3.3.3 Gráficos de variograma ajustado usando distancia y elevación.
png(file.path(ruta_datos,"figs/variograma_ajustado_distancia.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(v.dist.Jan.fit,plot.numbers=TRUE,cutoff=5000)
plot(v.dist.Apr.fit,plot.numbers=TRUE,cutoff=5000)
plot(v.dist.Jul.fit,plot.numbers=TRUE,cutoff=5000)
plot(v.dist.Oct.fit,plot.numbers=TRUE,cutoff=5000)
dev.off()## png
## 2
png(file.path(ruta_datos,"figs/variograma_ajustado_elevacion.png"), units="in", width=6, height=6, res=300)
par(mfrow=c(2,2))
plot(v.ele.Jan.fit,plot.numbers=TRUE,cutoff=5000)
plot(v.ele.Apr.fit,plot.numbers=TRUE,cutoff=5000)
plot(v.ele.Jul.fit,plot.numbers=TRUE,cutoff=5000)
plot(v.ele.Oct.fit,plot.numbers=TRUE,cutoff=5000)
dev.off()## png
## 2
4.3.3.4 Suma de cuadrado del error de variogramas.
rmse.v.Jan <- sqrt(sum(v.Jan $residuals)^2/nrow(v.Jan ))
rmse.v.Apr <- sqrt(sum(v.Apr $residuals)^2/nrow(v.Apr ))
rmse.v.Jul <- sqrt(sum(v.Jul $residuals)^2/nrow(v.Jul ))
rmse.v.Oct <- sqrt(sum(v.Oct $residuals)^2/nrow(v.Oct ))
rmse.v.dist.Jan <- sqrt(sum(v.dist.Jan $residuals)^2/nrow(v.dist.Jan ))
rmse.v.dist.Apr <- sqrt(sum(v.dist.Apr $residuals)^2/nrow(v.dist.Apr ))
rmse.v.dist.Jul <- sqrt(sum(v.dist.Jul $residuals)^2/nrow(v.dist.Jul ))
rmse.v.dist.Oct <- sqrt(sum(v.dist.Oct $residuals)^2/nrow(v.dist.Oct ))
rmse.v.ele.Jan <- sqrt(sum(v.ele.Jan $residuals)^2/nrow(v.ele.Jan ))
rmse.v.ele.Apr <- sqrt(sum(v.ele.Apr $residuals)^2/nrow(v.ele.Apr ))
rmse.v.ele.Jul <- sqrt(sum(v.ele.Jul $residuals)^2/nrow(v.ele.Jul ))
rmse.v.ele.Oct <- sqrt(sum(v.ele.Oct $residuals)^2/nrow(v.ele.Oct ))
rmse.v.dist.Jan.fit <- sqrt(sum(v.dist.Jan.fit $residuals)^2/nrow(v.dist.Jan.fit))
rmse.v.dist.Apr.fit <- sqrt(sum(v.dist.Apr.fit $residuals)^2/nrow(v.dist.Apr.fit))
rmse.v.dist.Jul.fit <- sqrt(sum(v.dist.Jul.fit $residuals)^2/nrow(v.dist.Jul.fit))
rmse.v.dist.Oct.fit <- sqrt(sum(v.dist.Oct.fit $residuals)^2/nrow(v.dist.Oct.fit))
rmse.v.ele.Jan.fit <- sqrt(sum(v.ele.Jan.fit $residuals)^2/nrow(v.ele.Jan.fit ))
rmse.v.ele.Apr.fit <- sqrt(sum(v.ele.Apr.fit $residuals)^2/nrow(v.ele.Apr.fit ))
rmse.v.ele.Jul.fit <- sqrt(sum(v.ele.Jul.fit $residuals)^2/nrow(v.ele.Jul.fit ))
rmse.v.ele.Oct.fit <- sqrt(sum(v.ele.Oct.fit $residuals)^2/nrow(v.ele.Oct.fit ))4.3.3.5 Interpretación y descripción de los variogramas.
4.4 Predicción espacial.
Cargar paquetes.
if(!require(gstat)) suppressMessages( suppressWarnings( install.packages("gstat") ))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))Cargar los datos.
archivos_rds <- list.files(file.path(ruta_datos,"rds"),
pattern="data_\\w+.rds")
archivos_rds <- archivos_rds [1:4]
archivos_rds## [1] "data_april.rds" "data_january.rds" "data_july.rds" "data_october.rds"data <- lapply(archivos_rds[1:4],function(file){
dat <- readRDS(file.path(ruta_datos,"rds",file))
id <- which(!is.na(dat$Tmax))
dat[id,]
})Cargamos de nuevo los datos mensuales para corregir datos faltantes.
dataJan <- readRDS(file.path(ruta_datos,"rds/data_january.rds"))
dataApr <- readRDS(file.path(ruta_datos,"rds/data_april.rds"))
dataJul <- readRDS(file.path(ruta_datos,"rds/data_july.rds"))
dataOct <- readRDS(file.path(ruta_datos,"rds/data_october.rds"))
dataJan <- dataJan [ !is.na(dataJan$dist), ]
dataApr <- dataApr [ !is.na(dataApr$dist), ]
dataJul <- dataJul [ !is.na(dataJul$dist), ]
dataOct <- dataOct [ !is.na(dataOct$dist), ]Modelos de regresión lineal.
mod <- lapply(data,function(dataMes){
lm(Tmax~TSS+dist+ele,dataMes)
})Vemos los valores de r-square, todos los r-square > 0,7.
lapply(mod,summary)## [[1]]
##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataMes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5353 -0.4794 0.0719 0.9757 2.5336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.546e+01 4.408e+00 3.507 0.00251 **
## TSS 3.553e-01 1.989e-01 1.786 0.09089 .
## dist 9.152e-06 1.498e-05 0.611 0.54891
## ele -3.206e-03 8.497e-04 -3.774 0.00139 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.651 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8568, Adjusted R-squared: 0.8329
## F-statistic: 35.89 on 3 and 18 DF, p-value: 8.35e-08
##
##
## [[2]]
##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataMes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9884 -0.5993 0.0650 0.9315 2.3007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.130e+01 3.761e+00 5.664 2.26e-05 ***
## TSS 2.495e-01 1.144e-01 2.181 0.04273 *
## dist 5.320e-05 1.427e-05 3.728 0.00154 **
## ele -6.103e-03 6.428e-04 -9.494 1.97e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.564 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8968, Adjusted R-squared: 0.8796
## F-statistic: 52.14 on 3 and 18 DF, p-value: 4.456e-09
##
##
## [[3]]
##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataMes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6094 -0.6500 -0.0589 0.5406 4.0711
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.068e+00 1.942e+00 4.153 0.000597 ***
## TSS 5.809e-01 1.359e-01 4.276 0.000455 ***
## dist -9.957e-06 1.457e-05 -0.683 0.503020
## ele 9.247e-04 9.855e-04 0.938 0.360528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.608 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.7781, Adjusted R-squared: 0.7411
## F-statistic: 21.04 on 3 and 18 DF, p-value: 4.111e-06
##
##
## [[4]]
##
## Call:
## lm(formula = Tmax ~ TSS + dist + ele, data = dataMes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.98716 -0.85230 0.09615 0.79237 1.59648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.711e+00 1.819e+00 5.338 4.50e-05 ***
## TSS 4.928e-01 7.762e-02 6.349 5.56e-06 ***
## dist -3.001e-06 1.135e-05 -0.264 0.7944
## ele -1.323e-03 7.294e-04 -1.813 0.0865 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.189 on 18 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9359, Adjusted R-squared: 0.9252
## F-statistic: 87.62 on 3 and 18 DF, p-value: 6.239e-11Graficamos Tmax vs los residuos.
lapply(mod,function(m){
plot(m$model$Tmax,m$residuals)
})
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL4.4.1 Modelos kriging para Julio.
4.4.1.1 Variogramas para Julio.
Ajustar modelos de variograma. Vamos a trabajar con el mes de Julio.
Extraigo del objeto lista 'data' el tercer elemento que corresponde al mes de Julio.
# dataJul <- data[[3]]
# dataJul <- st_transform(dataJul,32719)Variograma de Tmax.
v.jul <- variogram(Tmax~1,dataJul,
boundaries=c(0,5000,15000,seq(20000,150000,35000)))
plot(v.jul,plot.numbers=TRUE)
Variograma ajustado a Tmax.
v.jul.fit <- fit.variogram(v.jul,vgm(8,'Sph',50000))
plot(v.jul,v.jul.fit,plot.numbers=TRUE)
Variograma de los residuos.
v.jul.res <- variogram(Tmax~TSS+dist+ele,dataJul,
boundaries=c(0,seq(20000,40000,10000),
seq(60000,150000,20000)))
plot(v.jul.res,plot.numbers=TRUE)
Variograma ajustado a residuos.
v.jul.res.fit <- fit.variogram(v.jul.res,vgm(1,'Sph',60000))
plot(v.jul.res,v.jul.res.fit,plot.numbers=TRUE)
4.4.1.2 Variables explicatorias para Julio.
Cargamos la grilla con todos los predictores.
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))archivos_raster <- list.files(file.path(ruta_datos,"raster"),
pattern="dist_costa_biobio.tif|elevacion_biobio.tif|MOD11C3_A2017_tss_stack.tif")
archivos_raster## [1] "dist_costa_biobio.tif" "elevacion_biobio.tif" "MOD11C3_A2017_tss_stack.tif"Cargamos raster de TSS para el mes de Julio.
TSS <- raster(file.path(ruta_datos,"raster",archivos_raster[3]),3)Juntamos en un stack los raster de TSS, dist y ele.
grid <- stack(file.path(ruta_datos,"raster",archivos_raster[c(1,2)]))Cortamos para igualar los extent.
TSS <- crop(TSS,grid)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
TSS <- TSS-273.15Creamos el stack con los tres predictores.
grid <- stack(TSS,grid)
names(grid) <- c('TSS','dist','ele')
plot(grid)
grid.ok <- as(grid,'SpatialGridDataFrame')Transformamos coordenadas de raster a 32719.
# grid <- projectRaster(grid,crs=CRS('+init=epsg:32719'))# crs(grid) <- NA
# st_crs(dataJul) <- NAInverso de la distancia
res.idw <- krige(Tmax~1,as(dataJul,'Spatial'),grid.ok)## [inverse distance weighted interpolation]res.idw <- stack(res.idw)
plot(res.idw)
4.4.1.3 Kriging Ordinario para Julio.
res.ok <- krige(Tmax~1,as(dataJul,'Spatial'),grid.ok,model=v.jul.fit)## [using ordinary kriging]res.ok <- stack(res.ok)plot(res.ok)
4.4.1.4 Regresión kriging para Julio.
res.rk <- krige(Tmax~TSS+dist+ele,as(dataJul,'Spatial'),
grid.ok,model=v.jul.res.fit)## [using universal kriging]res.rk <- stack(res.rk)
plot(res.rk)
Visualizamos usando tmap.
if(!require(tmap)) suppressMessages( suppressWarnings( install.packages("tmap") ))
# tmap_mode('view')Inverso de la distancia asigno CRS a res.idw
# crs(res.idw) <- CRS('+init=epsg:32719')
# st_crs(dataJul) <- 32719tm_shape(res.idw[[1]]) +
tm_raster() +
tm_shape(dataJul) +
tm_dots()
Kriging ordinario asigno CRS a res.idw
# crs(res.ok) <- CRS('+init=epsg:32719')tm_shape(res.ok[[1]]) +
tm_raster() +
tm_shape(dataJul) +
tm_dots()
4.4.1.5 Ordinary Kriging con validación cruzada para Julio.
res.ok.cv <- krige.cv(Tmax~1,as(dataJul,'Spatial'),
grid.ok,model=v.jul.fit,nfold=10)r2 <- 1- sum(res.ok.cv$residual^2)/sum((res.ok.cv$observed-mean(res.ok.cv$observed))^2)4.4.1.6 Regresión kriging con validación cruzada para Julio.
res.rk.cv <- krige.cv(Tmax~TSS+dist+ele,as(dataJul,'Spatial'),
grid.ok,model=v.jul.res.fit,nfold=10)r2 <- 1- sum(res.rk.cv$residual^2)/sum((res.rk.cv$observed-mean(res.rk.cv$observed))^2)4.4.2 Modelos kriging para Octubre.
4.4.2.1 Variogramas para Octubre.
Ajustar modelos de variograma. Vamos a trabajar con el mes de Octubre.
Extraigo del objeto lista 'data' el tercer elemento que corresponde al mes de Octubre.
# dataOct <- data[[4]]
# dataOct <- st_transform(dataOct,32719)Variograma de Tmax.
v.oct <- variogram(Tmax~1,dataOct,
boundaries=c(0,10000,20000,seq(35000,150000,20000)))
plot(v.oct,plot.numbers=TRUE)
Variograma ajustado a Tmax.
v.oct.fit <- fit.variogram(v.oct,vgm(15,'Exp',120000))
plot(v.oct,v.oct.fit,plot.numbers=TRUE)
Variograma de los residuos.
v.oct.res <- variogram(Tmax~TSS+dist+ele,dataOct,
boundaries=c(0,5000,10000,
seq(15000,150000,20000)))
plot(v.oct.res,plot.numbers=TRUE)
Variograma ajustado a residuos.
v.oct.res.fit <- fit.variogram(v.oct.res,vgm(1,'Exp',60000))
plot(v.oct.res,v.oct.res.fit,plot.numbers=TRUE)
4.4.2.2 Variables explicatorias para Octubre.
Cargamos la grilla con todos los predictores.
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))archivos_raster <- list.files(file.path(ruta_datos,"raster"),
pattern="dist_costa_biobio.tif|elevacion_biobio.tif|MOD11C3_A2017_tss_stack.tif")
archivos_raster## [1] "dist_costa_biobio.tif" "elevacion_biobio.tif" "MOD11C3_A2017_tss_stack.tif"Cargamos raster de TSS para el mes de Octubre.
TSS <- raster(file.path(ruta_datos,"raster",archivos_raster[3]),4)Juntamos en un stack los raster de TSS, dist y ele.
grid <- stack(file.path(ruta_datos,"raster",archivos_raster[c(1,2)]))Cortamos para igualar los extent.
TSS <- crop(TSS,grid)Aplicamos resta de 273.15 como conversión de grados Kelvin a Celsius.
TSS <- TSS-273.15Creamos el stack con los tres predictores.
grid <- stack(TSS,grid)
names(grid) <- c('TSS','dist','ele')
plot(grid)
grid.ok <- as(grid,'SpatialGridDataFrame')Transformamos coordenadas de raster a 32719.
# grid <- projectRaster(grid,crs=CRS('+init=epsg:32719'))# crs(grid) <- NA
# st_crs(dataOct) <- NAInverso de la distancia.
res.idw <- krige(Tmax~1,as(dataOct,'Spatial'),grid.ok)## [inverse distance weighted interpolation]res.idw <- stack(res.idw)
plot(res.idw)
4.4.2.3 Kriging Ordinario para Octubre.
res.ok <- krige(Tmax~1,as(dataOct,'Spatial'),grid.ok,model=v.oct.fit)## [using ordinary kriging]res.ok <- stack(res.ok)plot(res.ok)
4.4.2.4 Regresión kriging para Octubre.
res.rk <- krige(Tmax~TSS+dist+ele,as(dataOct,'Spatial'),
grid.ok,model=v.oct.res.fit)## [using universal kriging]res.rk <- stack(res.rk)
plot(res.rk)
Visualizamos usando tmap.
if(!require(tmap)) suppressMessages( suppressWarnings( install.packages("tmap") ))
# tmap_mode('view')Inverso de la distancia asigno CRS a res.idw
# crs(res.idw) <- CRS('+init=epsg:32719')
# st_crs(dataOct) <- 32719tm_shape(res.idw[[1]]) +
tm_raster() +
tm_shape(dataOct) +
tm_dots()
Kriging ordinario asigno CRS a res.idw
# crs(res.ok) <- CRS('+init=epsg:32719')tm_shape(res.ok[[1]]) +
tm_raster() +
tm_shape(dataOct) +
tm_dots()
4.4.2.5 Ordinary kriging con validación cruzada para Octubre.
res.ok.cv <- krige.cv(Tmax~1,as(dataOct,'Spatial'),
grid.ok,model=v.oct.fit,nfold=10)r2 <- 1- sum(res.ok.cv$residual^2)/sum((res.ok.cv$observed-mean(res.ok.cv$observed))^2)4.4.2.6 Regresión kriging con validación cruzada para Octubre.
res.rk.cv <- krige.cv(Tmax~TSS+dist+ele,as(dataOct,'Spatial'),
grid.ok,model=v.oct.res.fit,nfold=10)r2 <- 1- sum(res.rk.cv$residual^2)/sum((res.rk.cv$observed-mean(res.rk.cv$observed))^2)4.5 Evaluación de la predicción espacial.
4.5.0.1 Regresión alternativa para replicar.
Alternativa de prueba: ajustar el variograma con la Tmax estimada obtenida del predictor de TSS y luego utilizarlo para interpolar la Tmax observada.
TSS <- grid[[1]]i <- sample(ncell(TSS),ncell(TSS)*.2)
TSS.ex <- TSS[i]
coords <- xyFromCell(TSS,i)dataTSS <- data.frame(coords,TSS.ex)
id <- which(!is.na(dataTSS[,3]))
dataTSS <- dataTSS[id,]dataTSS <- st_as_sf(dataTSS,coords = c('x','y'))Cambiamos nombre de la variable para que R crea que es Tmax, pero en realidad es TSS.
names(dataTSS)[1] <- 'Tmax'
var.tss <- variogram(Tmax~1,dataTSS)
plot(var.tss)
var.tss.fit <- fit.variogram(var.tss,vgm(30,'Sph',150000))
plot(var.tss,var.tss.fit,plot.numbers = TRUE)
Interpolar con OK usando el variograma del muestreo de TSS.
grid.ok <- as(grid,'SpatialGridDataFrame')# crs(grid.ok) <- NA
# st_crs(dataOct) <- NAres.ok <- krige(Tmax~1,as(dataOct,'Spatial'),grid.ok,model=var.tss.fit)## [using ordinary kriging]res.ok <- stack(res.ok)plot(res.ok)
4.5.0.2 Ordinary kriging + cross validation para Octubre.
res.ok.cv <- krige.cv(Tmax~1,as(dataOct,'Spatial'),
grid.ok,model=var.tss.fit,nfold=10)r2 <- 1- sum(res.ok.cv$residual^2)/sum((res.ok.cv$observed-mean(res.ok.cv$observed))^2)4.5.0.3 Regresión kriging + cross validation para Octubre.
res.rk.cv <- krige.cv(Tmax~TSS+dist+ele,as(dataOct,'Spatial'),
grid.ok,model=vgm(1.73,'Sph',146065),nfold=10)r2 <- 1- sum(res.rk.cv$residual^2)/sum((res.rk.cv$observed-mean(res.rk.cv$observed))^2)rmse <- sqrt(sum(res.rk.cv$residual^2)/36)4.5.0.4 Predicción con regresion-kriging para Octubre.
res.rk <- krige(Tmax~TSS+dist+ele,as(dataOct,'Spatial'),
grid.ok,model=vgm(1.73,'Sph',146065))## [using universal kriging]res.rk <- stack(res.rk)
plot(res.rk)
Visualizamos usando tmap.
if(!require(tmap)) suppressMessages( suppressWarnings( install.packages("tmap") ))
# tmap_mode('view')Inverso de la distancia. Asignamos CRS a res.idw.
# crs(res.rk) <- CRS('+init=epsg:32719')
# st_crs(dataOct) <- 32719tm_shape(res.rk[[2]]) +
tm_raster(palette='YlOrRd') +
tm_shape(dataOct) +
tm_dots()
5 Discusión.
En términos de la regresión lineal observamos que ésta ignora fenómenos de autocorrelación presente y que al momento de usar una herramienta como el variograma uno puede encontrar una gran variedad de opciones para poder identificar y modelar dicho fenómeno y de este modo conocer con mayor precisión la muestra que se está analizando.
Hemos observado que la regresión kriging entregó en algunos meses una curva con un punto de inflección más agudo, lo que nos indica que en virtud de esos parámetros, se pueden realizar pronósticos con una base más segura.
6 Conclusión.
Hemos aprendido con las técnicas de geoestadística que para poder cuantificar entidades geográficas es necesario abordar aspectos tales como dependencia espacial. Hemos podido hacer una aproximación a poder diagnosticar y concer algunas estructuras de correlación espacial que puede presentarse en fenómenos conocidos de la naturaleza. Por su parte la regresión kriging nos brinda una mejor aproximación para poder realizar predicción espacial de lo que es un modelo de regresión lineal ya que permite abordar esquemas de dependencia o correlación presentes en variables de interés meteorológicos o similar.