Actividad 5
Nota. Este documento se puede ver mejor en html, puesto que contiene mapas dinámicos en https://rpubs.com/mosorio/982481
Ejemplo aplicado de Muestreo espacial
La geoestadística son las herramientas y técnicas que sirven para analizar y predecir los valores de una variable o más que se muestra distribuida en el espacio o en el tiempo de una forma continua.
En el muestreo espacial, se recolectan muestras en un marco bidimensional (latitud y longitud). El esquema de muestreo está diseñado para maximizar la probabilidad de incluir la variación espacial de la población de estudio. Existen diferentes métodos para recolectar una muestra espacial, en esta actividad consideraremos tres métodos de muestreo ya revisados en actividades anteriores:
Muestreo aleatorio simple: selección aleatoria de puntos de muestra en una región de estudio, donde cada ubicación tiene la misma oportunidad de ser muestreada.
Muestreo sistemático: el muestreo sistemático es un tipo de método de muestreo probabilístico en el que los miembros de la muestra de una población más grande se seleccionan de acuerdo con un punto de partida aleatorio pero con un intervalo periódico fijo.
Muestreo estratificado: la región de estudio se divide en grupos, también conocidos como estratos, según una característica colectiva de la región de estudio (elevación, uso del suelo, etc.). Luego, para cada estrato, se recolecta una muestra aleatoria espacial y se combina en una muestra.
De los primeros paquetes en R que se desarrollaron para proveer clases y métodos para datos espaciales fueron sp y rgdal. En un primer momento los usaremos para graficar un área y seleccionar puntos aleatorios de esta área.
Lo siguiente es un ejemplo de cómo generar sitios de muestreo usando el paquete sp. El mapa que se usará es la del estado de Oregon en Estados Unidos.
# paquetes usados
library(sp)
library(rgdal)
# leer el archivo shape de Oregon
x <- readOGR('orotl.shp', layer='orotl')## OGR data source with driver: ESRI Shapefile
## Source: "/home/mucio/Master Estadística/Encuestas por muestreo. Aplicaciones económicas, sociales y medioambientales/Tema 5. Metodos de muestreo para encuestas sociales y economicas/espacial/orotl.shp", layer: "orotl"
## with 36 features
## It has 1 fields# graficar el mapa de Oregon
plot(x, main = "Estado de Oregon")
La función spsample del paquete sp genera ubicaciones de puntos dentro de un área cuadrada, una cuadrícula, un polígono o en una línea espacial, utilizando métodos de muestreo regulares o aleatorios; los métodos utilizados asumen que la geometría utilizada no es esférica, por lo que los objetos deben estar en coordenadas planas.
Muestreo simple aleatorio
Las siguientes líneas generan 100 coordenadas espaciales aleatorias simples en el estado de Oregon.
# ajustar la semilla aleatoria
set.seed(1234)
# generar puntos aleatorios
puntos <- spsample(x, n=100, type='random')
# ver las latitudes y longitudes elegidas
head(puntos@coords,20)## x y
## [1,] -121.5679 42.59138
## [2,] -123.6247 43.85583
## [3,] -118.5109 44.42466
## [4,] -121.7706 45.67125
## [5,] -122.9818 44.64392
## [6,] -118.1121 42.90065
## [7,] -122.4139 45.15884
## [8,] -121.5173 42.38155
## [9,] -123.0036 45.28447
## [10,] -118.9261 42.43491
## [11,] -123.5828 43.54443
## [12,] -121.2642 44.66619
## [13,] -122.9676 43.06988
## [14,] -122.5919 43.42715
## [15,] -121.6687 42.78947
## [16,] -117.9155 45.41972
## [17,] -118.4274 44.80661
## [18,] -119.1585 44.79324
## [19,] -120.6119 42.53117
## [20,] -123.0164 45.73668# graficar los puntos en el estado
plot(x)
points(puntos, col = 'red', pch = 3, cex = 0.5)
Muestreo espacial sistemático
El segundo método de muestreo que usaremos se llama muestreo sistemático. Para hacerlo, primero se crea una cuadrícula regular en la parte superior del área de estudio y luego seleccionar el punto central de cada celda de la cuadrícula. Una vez que se tienen todos los puntos centrales, esto será el marco de muestreo. Entonces se puede seleccionar del marco de muestreo basado en un método de selección sistemático (como cualquier otro punto y así sucesivamente).
# ajustar la semilla aleatoria
set.seed(1234)
# generar puntos aleatorios
puntos_sist <- spsample(x, n=100, type='regular')
# ver las latitudes y longitudes elegidas
head(puntos_sist@coords,20)## x1 x2
## [1,] -123.9673 42.31807
## [2,] -123.4366 42.31807
## [3,] -122.9058 42.31807
## [4,] -122.3751 42.31807
## [5,] -121.8443 42.31807
## [6,] -121.3136 42.31807
## [7,] -120.7829 42.31807
## [8,] -120.2521 42.31807
## [9,] -119.7214 42.31807
## [10,] -119.1906 42.31807
## [11,] -118.6599 42.31807
## [12,] -118.1292 42.31807
## [13,] -117.5984 42.31807
## [14,] -117.0677 42.31807
## [15,] -124.4980 42.84881
## [16,] -123.9673 42.84881
## [17,] -123.4366 42.84881
## [18,] -122.9058 42.84881
## [19,] -122.3751 42.84881
## [20,] -121.8443 42.84881# graficar los puntos en el estado
plot(x)
points(puntos_sist, col = 'red', pch = 3, cex = 0.5)
Muestreo estratificado
La función cuando se requiere el muestreo aleatorio estratificado, toma una muestra aleatoria uniforme de tamaño k de cada uno de los m estratos o subáreas, por lo que n = km.
set.seed(1234)
plot(x)
points(spsample(x, n=100, type='stratified'), col='red', pch=3, cex=0.5)
Ahora veamos como funcionan otros paquetes para la selección espacial de sitios. Se usaran datos de biomasa sobre el suelo. La biomasa aérea (AGB) se define como “la masa seca sobre el suelo de materia viva o muerta de formas de vida de árboles o arbustos (leñosas), expresada como masa por unidad de área”.
Se aplicaran diferentes métodos de muestreo para estimar la biomasa aérea en la isla de Vancouver (área de estudio) y se comparará la media muestral estimada con tres métodos de muestreo (simple aleatorio, sistemático y estratificdo) con la media real de la población.
Ahora se usarán otros paquetes más:
library(raster)
library(rgdal)
library(sf)## Linking to GEOS 3.11.1, GDAL 3.6.0, PROJ 9.1.0; sf_use_s2() is TRUElibrary(sp)
library(leaflet)
library(tidyverse)## ── Attaching packages
## ────────────────────────────────── tidyverse
## 1.3.2.9000 ──## ✔ ggplot2 3.4.0 ✔ dplyr 1.0.10
## ✔ tibble 3.1.8 ✔ stringr 1.4.1
## ✔ tidyr 1.2.1 ✔ forcats 0.5.2
## ✔ readr 2.1.3 ✔ lubridate 1.9.0
## ✔ purrr 0.3.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()library(terra)## terra 1.6.47
##
## Attaching package: 'terra'
##
## The following object is masked from 'package:tidyr':
##
## extract
##
## The following object is masked from 'package:rgdal':
##
## projectGeneralmente no se tiene acceso a la población. Por eso se usan métodos de muestreo. Sin embargo, en esta actividad tenemos una capa de mapa que incluye datos AGB. Supondremos que esta capa del mapa es nuestra población. Entonces, si se calcula el valor promedio de los datos de esta capa, se tiene la media de la población \(\mu\).
Entonces, si seleccionamos un conjunto de datos de esta capa del mapa AGB, podemos asumirlo como una muestra recopilada de la población.
Este conjunto de datos en particular se descarga desde (https://open.canada.ca/data/en/dataset/ec9e2659-1c29-4ddb-87a2-6aced147a990) y puede leer más sobre él en el enlace proporcionado. Estos datos son de 2001 a 2011 y están disponibles para todo Canadá.
El formato de los datos es raster y usaremos para graficarlos al paquete leaflet
# lectura de datos raster
biomass <- raster("https://www.dropbox.com/s/zvzqb4a0qvhulye/bc_NFI_MODIS250m_2011_kNN_Structure_Biomass_TotalLiveAboveGround_v1_victoria_4326.tif?dl=1")
#la paleta de colores
biomass_pallet <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(biomass),na.color = "transparent")
# gráficar en mapa interactivo
m <- leaflet() %>%
addTiles() %>%
addRasterImage(biomass, opacity = 0.8) %>%
addLegend(pal = biomass_pallet, values = values(biomass),title = "Biomasa total\n Live Above Ground")
mLos valores de Biomasa (en total 253948) se pueden ver usando:
# Se remueven los valores NA
p_data <- values(biomass)[!is.na(values(biomass))]
head(p_data) ## [1] 12.3500 207.4633 205.9867 207.4633 209.9700 209.9700length(p_data)## [1] 253948La distribución de frecuencias de la biomasa no es distribuida de manera normal.
hist(p_data, main = "Biomasa", xlab = "t/ha")
Si consideramos a estos datos como la población, entonces los parámetros de localización central y de variación son:
p_mean <- mean(p_data)
p_mean## [1] 229.6976p_sd <- sqrt(sum((p_data - p_mean)^2) / (length(p_data)))
p_sd## [1] 120.597Ahora, veamos las estimaciones de estos parámetros con los tres tipos de muestreo arriba mencionados.
Muestreo aleatorio simple
El paquete raster tiene también funciones para la selección de muestras. Considere que se seleccionarán de manera aleatoria 100 sitios de la región.
La siguiente instrucción genera 100 sitios aleatorios del marco de muestreo que ya tenemos.
# ajustamos la semilla
set.seed(1234)
# la muestra de tamaño 100
muestra_1 <- sampleRandom(biomass, size=100, cells=TRUE, sp=TRUE)
head(muestra_1@data)# las ubicaciones de estos 100 sitios
head(muestra_1@coords)## x y
## [1,] -126.5702 50.43309
## [2,] -125.0068 49.39217
## [3,] -125.0495 49.42279
## [4,] -124.9396 49.38962
## [5,] -124.5060 48.83090
## [6,] -125.6297 49.34625En el mapa, los 100 sitios
mapa_muestra_1 <- m %>%
leaflet::addCircles(muestra_1@coords[,1],muestra_1@coords[,2],popup = paste("Pixel Value:",muestra_1$bc_NFI_MODIS250m_2011_kNN_Structure_Biomass_TotalLiveAboveGround_v1_victoria_4326.tif.dl.1),color="red")
mapa_muestra_1El cálculo del promedio y desviación estándar muestral
muestra_1_media <- mean(muestra_1$bc_NFI_MODIS250m_2011_kNN_Structure_Biomass_TotalLiveAboveGround_v1_victoria_4326)
muestra_1_media## [1] 222.5427muestra_1_sd <- sd(muestra_1$bc_NFI_MODIS250m_2011_kNN_Structure_Biomass_TotalLiveAboveGround_v1_victoria_4326)
muestra_1_sd## [1] 110.6047El error de muestreo en la estimación del promedio a través del muestreo simple aleatorio
e1 <- abs(p_mean-muestra_1_media)
e1## [1] 7.154891Muestreo aleatorio sistemático
Como se había dicho en este muestreo se genera una cuadricula y luego se selecciona el centro de esta. Para generar la cuadricula regular se usa la función st_make_grid. La primera entrada de esta función es el archivo ráster sobre el que desea dibujar una cuadrícula. El segundo parámetro es el tipo de cuadrícula de salida y el tercer parámetro de esta función es el número de filas y columnas de la cuadrícula. En este caso vamos a tener una cuadrícula de 20 por 20 (para apreciar los puntos mejor).
# cuadricula de 20 * 20 sobre el área de estudio
grid <- st_make_grid(biomass,what = "polygons",n = c(20, 20))
map_grid <- m %>% addPolygons(data = grid)
map_gridAhora, se generan los centros de estas cuadriculas
# la misma función pero ahora con centros como argumento
grid <- st_make_grid(biomass,what = "centers",n = c(20, 20))
map_grid <- m %>%
addCircles(data = grid)
map_gridComo se puede ver, la cuadrícula de muestreo tiene celdas que están fuera del área de estudio. Para limpiar esas celdas. Además se debe asignar un valor de ráster debajo de cada punto a su punto representativo.
grid_sf <- st_as_sf(grid)
grid_points_with_values <- extract(biomass,grid_sf)
grid_sf <- st_sf(value=grid_points_with_values,grid)
grid_sf <- grid_sf[!is.na(grid_sf$value),] #remover NA
map_sample_2 <- m %>% addCircles(data = grid_sf,popup = paste("Biomas Value:", grid_sf$value),color="yellow")
map_sample_2Para seleccionar 50 puntos de éstos primero se calcula la constante de salto sistematico
# el total de puntos de la cuadrícula entre el tamaño de muestra 50
sample_2_interval <- round(length(grid_sf$value)/50)
sample_2_interval## [1] 2# select sample items based on the sample interval we just calculated
grid_sf_grid <- grid_sf$grid[1:(length(grid_sf$grid)/sample_2_interval)*sample_2_interval]
grid_sf_value <- grid_sf$value[1:(length(grid_sf$value)/sample_2_interval)*sample_2_interval]
#MOstrar puntos en mapa
map_sample_2 %>% addCircles(data=grid_sf_grid,color="red")El estimador del promedio muestral con el muestreo sistemático
muestra_2_media <- mean(grid_sf_value)
muestra_2_media## [1] 232.4691muestra_2_sd <- sd(grid_sf_value)
muestra_2_sd## [1] 100.4802El error de estimación
e2 <- abs(p_mean-muestra_2_media)
e2## [1] 2.771498Muestreo estratificado
Para estratificar se usará la elevación de los sitios. Existe información sobre la elevación de los sitios y esta información se almacena en archivos DEM.
dem <- raster("https://www.dropbox.com/s/nrxzy6bemzg4tks/dem_lowresolution.tif?dl=1")
# se reemplazan los valores con elevaciones menores o iguales que cero con NA
values(dem)[values(dem) <= 0] = NA
# graficar en el mapa
dem_pallet <- colorNumeric(c("#008435", "#33cc00", "#f4f071","#f4bd45","#99642b","#ffffff"), values(dem),na.color = "transparent")
map_dem <- leaflet() %>%
addTiles() %>%
addRasterImage(dem, colors = dem_pallet, opacity = 0.8) %>%
addLegend(pal = dem_pallet, values = values(dem),title = "DEM (elevación en m)") ## Warning in colors(.): Some values were outside the color scale and will be
## treated as NAmap_demAhora, si se definen los estratos
- Estrato 1: Elevaciones menores o iguales que 900 m
- Estrato 2: Elevaciones mayores que 900 m
estrato_1 <- dem
# remplazar valores <= 900 con 1 y con NA a los > 900
values(estrato_1)[values(estrato_1) <= 900] = 1
values(estrato_1)[values(estrato_1) > 900] = NA
estrato_2 <- dem
# remplazar valores > 900 con 2 y otros con NA
values(estrato_2)[values(estrato_2) <= 900] = NA
values(estrato_2)[values(estrato_2) > 900] = 2 En el mapa
strata_bins <- colorBin(c("#008435", "#99642b"),bins = 2, c(1,2),na.color = "transparent")
leaflet() %>%
addTiles() %>%
addRasterImage(estrato_1, colors = strata_bins, opacity = 0.8)%>%
addRasterImage(estrato_2, colors = strata_bins, opacity = 0.8)%>%
addLegend(pal = strata_bins, values = c(1,2),title = "Estratos 1 y 2") ## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NAAhora se realiza el muestreo aleatorio en cada estrato (50 en cada estrato)
# limpieza de NA
values(estrato_1)[is.na(values(biomass))] = NA
values(estrato_2)[is.na(values(biomass))] = NA
# muestreo aleatorio en cada estrato
muestra_3_est1 <- sampleRandom(estrato_1, size=round(50), cells=TRUE, sp=TRUE)
muestra_3_est2 <- sampleRandom(estrato_2, size=round(50), cells=TRUE, sp=TRUE)
# extraer los valores de cada estrato
values_st1 <- values(biomass)[muestra_3_est1$cell]
values_st2 <- values(biomass)[muestra_3_est2$cell]
# en el mapa
map_sample_3 <- m %>%
leaflet::addCircles(muestra_3_est1@coords[,1],muestra_3_est1@coords[,2],popup = paste("Pixel Value (Strata 1):",values_st1),color = "red") %>%
leaflet::addCircles(muestra_3_est2@coords[,1],muestra_3_est2@coords[,2],popup = paste("Pixel Value (Strata 2):",values_st2),color = "orange")
map_sample_3La estimación del promedio con el muestreo estratificado
# media y desviación estándar (50% en cda estrato de los valores)
muestra_3_media <- 1/2*mean(values_st1) + 1/2*mean(values_st2)
muestra_3_media## [1] 209.0838El error de estimación
e3 <- abs(p_mean-muestra_3_media)
e3## [1] 20.61382Con el muestreo simple aleatorio el error es 7.1548914, con el muestreo sistemático 2.7714977 y con el muestreo estratificado 20.613824.
Documentos tomados como referencia:
https://rpubs.com/majidhojati/lab4
https://rstudio-pubs-static.s3.amazonaws.com/200263_e77d00d6b6b24aa8890c8c4f074bcdff.html