para - poligono en r

Trazando datos interpolados en el mapa

Question

Tengo datos de encuestas sobre la riqueza de especies que se tomaron en varios sitios en la bahía de Chesapeake, Estados Unidos, y me gustaría presentar gráficamente los datos como un "mapa de calor".

Tengo un marco de datos de coordenadas lat / long y valores de riqueza, que convertí en un SpatialPointsDataFrame y usé la función autoKrige() del paquete automap para generar los valores interpolados.

Primero, ¿alguien puede comentar si estoy implementando correctamente la función autoKrige() ?

En segundo lugar, tengo problemas para trazar los datos y superponer un mapa de la región. Alternativamente, ¿podría especificar la cuadrícula de interpolación para reflejar los bordes de la bahía (como se sugiere here )? ¿Alguna idea sobre cómo podría hacer eso y dónde podría obtener esa información? Suministrar la grilla a autoKrige() parece bastante fácil.

EDITAR: ¡Gracias a Paul por su súper útil publicación! Esto es lo que tengo ahora. Tiene problemas para hacer que ggplot acepte tanto los datos interpolados como la proyección del mapa:

require(rgdal)
require(automap)
#Generate lat/long coordinates and richness data
set.seed(6)
df=data.frame(
lat=sample(seq(36.9,39.3,by=0.01),100,rep=T),
long=sample(seq(-76.5,-76,by=0.01),100,rep=T),
fd=runif(10,0,10))
initial.df=df
#Convert dataframe into SpatialPointsDataFrame
coordinates(df)=~long+lat
#Project latlong coordinates onto an ellipse
proj4string(df)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
#+proj = the type of projection (lat/long)
#+ellps and +datum = the irregularity in the ellipse represented by planet earth
#Transform the projection into Euclidean distances
project_df=spTransform(df, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")) #projInfo(type="proj")
#Perform the interpolation using kriging
kr=autoKrige(fd~1,project_df)
#Extract the output and convert to dataframe for easy plotting with ggplot2
kr.output=as.data.frame(kr$krige_output)
#Plot the output
#Load the map data for the Chesapeake Bay
cb=data.frame(map("state",xlim=range(initial.df$long),ylim=range(initial.df$lat),plot=F)[c("x","y")])
ggplot()+
geom_tile(data=kr.output,aes(x=x1,y=x2,fill=var1.pred))+
geom_path(data=cb,aes(x=x,y=y))+
coord_map(projection="mercator")

Answer 1

Tengo una serie de comentarios en tu publicación:

Usando kriging

Veo que está utilizando geoestadística para construir su mapa de calor. También podría considerar otras técnicas de interpolación, como splines (p. Ej., Splines de placa delgada en el paquete de campos). Estos hacen menos suposiciones sobre los datos (por ejemplo, estacionariedad), y también pueden visualizar sus datos muy bien. La reducción en el número de suposiciones podría ayudar en caso de que lo envíe a un diario, entonces tendrá menos para explicar a los revisores. También puede comparar algunas técnicas de interpolación si lo desea, consulte un informe que escribí para obtener algunos consejos.

Proyección de datos

Veo que está utilizando coordenadas de latitud larga para kriging. Edzer Pebesma (autor de gstat ) comentó que no hay modelos de variograma que sean adecuados para las coordenadas lat lon. Esto se debe a que en latón las distancias no son rectas (es decir, Euclidean ), sino en una esfera (es decir, distancias de gran círculo ). No hay funciones de covarianza (o modelos de variograma) que sean válidas para coordenadas esféricas. Recomiendo proyectarlos utilizando spTransform desde el paquete rgdal antes de usar automap.

El paquete rgdal usa la biblioteca de proyección proj4 para realizar los cálculos. Para proyectar sus datos, primero debe definir su proyección:

proj4string(df) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

La cadena proj4 en el lado derecho de la expresión anterior define el tipo de proyección ( +proj ), las elipses que se usaron ( +ellps ) y el dato ( +datum ). Para entender lo que significan estos términos, debes imaginar la Tierra como una papa. La Tierra no es perfectamente esférica, esto está definido por las elipses. Tampoco la Tierra es un elipsoide perfecto, pero la superficie es más irregular. Esta irregularidad está definida por el dato. Ver también este artículo en Wikipedia .

Una vez que tenga la proyección definida, puede usar spTransform :

project_df = spTransform(df, CRS("+proj= etcetc"))

donde CRS ("+ proj, etc.") define la proyección del objetivo. La proyección adecuada depende de su ubicación geográfica y del tamaño de su área de estudio.

Trazando con ggplot2

Para agregar polígonos o polilíneas a ggplot, consulte la documentación de coord_map . Esto incluye un ejemplo de uso del paquete de maps para trazar límites de país. Si necesita cargar, por ejemplo, shapefiles para su área de estudio, puede hacerlo usando rgdal . Recuerde que ggplot2 funciona con data.frame''s, no SpatialPolygons . Puede transformar SpatialPolygons en data.frame usando:

poly_df = fortify(poly_Spatial)

Ver también esta función que creé para trazar grillas espaciales. Funciona directamente en SpatialGrids / Pixels. Tenga en cuenta que necesita obtener uno o dos archivos adicionales de ese repositorio ( continuousToDiscrete ).

Crear una grilla de interpolación

Creé automap para generar una grilla de salida cuando no se especificó ninguna. Esto se hace creando un casco convexo alrededor de los puntos de datos y muestreando 5000 puntos dentro de él. Los límites del área de predicción y el número de puntos muestreados en ella (y, por lo tanto, la resolución) es bastante arbitraria. Para una aplicación específica, la forma del área de predicción puede derivarse de un polígono, usando spsample para muestrear puntos dentro del polígono. Cuantos puntos analizar, y por lo tanto la resolución, depende de dos cosas:

• el tipo de datos que tiene, por ejemplo, si sus datos son muy sencillos, no tiene mucho sentido elevar la resolución realmente alta en comparación con esta suavidad. Alternativamente, si sus datos tienen muchas líneas a pequeña escala, necesita una resolución alta. Esto solo es posible si tienes las observaciones para apoyar esta alta resolución.
• la densidad de datos. Si sus datos son más densos, puede elevar la resolución.

Si usa su mapa interpolado para análisis posteriores, es importante obtener la resolución correcta. Si usas el mapa solo para fines de visuatización, esto es menos importante. Sin embargo, tenga en cuenta que, en ambos casos, una resolución demasiado alta puede ser engañosa en cuanto a la precisión de sus predicciones, y que una resolución demasiado baja no hace justicia a los datos.