Guardar/exportar niveles de contorno específicos para un KDE (paquete 'ks' en R) como shapefile de ESRI

Question

Estoy realizando un análisis del núcleo (KDE) de un individuo con el paquete 'ks' en R. Me gustaría evitar el uso de Geospatial Modelling Environment (GME) para implementar la salida en ArcGIS, porque estoy usando QGIS que no es compatible con GME. Me gustaría ejecutar el análisis en R (véase el ejemplo siguiente) e importar el resultado como un shapefile en QGIS.

gps <- read_csv("data/gps.csv")
gps<-data.frame(cbind(gps$x,gps$y))
colnames(gps)<-c('x','y')
gps_hpi<-Hpi(x=gps,pilot="samse",pre="scale")
gps_kde<-kde(x=gps,h=gps_hpi,compute.cont=TRUE)
contourLevels(gps_kde,cont=c(50,75,95))
contourSizes(gps_kde,cont=c(50,75,95),approx=TRUE)
plot(gps_kde, display="filled.contour", cont=c(50,75,95))

Estos contornos son exactamente lo que quiero exportar como un shapfile, pero todavía no hay ninguna posibilidad razonable descrita en internet para hacerlo.

La adaptación de la fórmula de @Jeffrey Evans dio como resultado un KDE-plot que difiere significativamente de la función original plot.kde.

r<-raster(extent(440000,465000,5655000,5680000),nrows=nrow(gps_kde$estimate), ncols=ncol(gps_kde$estimate))  
r[]<-gps_kde$estimate
r<-flip(r,direction='y')
rcont<-rasterToContour(r,levels=gps_kde$cont)
(cont.values<-gps_kde$cont[grep(paste(c("50","75","95"),collapse="|"),names(gps_kde$cont))])
rcont.gt50<-rcont[which(rcont@data[,1] %in% cont.values),]
par(mfrow=c(1,1))
plot(gps_kde,display="image", main="KDE")
plot(r,main="raster of KDE")
plot(rcont,main="All contours (1-99%)")
plot(rcont.gt50, main="50%, 75% and 95% volume contours")

Con este resultado no es posible realizar más análisis. También el método de @Juan Antonio Roldán Díaz dio resultados no aplicables.

hts<-contourLevels(gps_kde, prob = c(0.5, 0.75, 0.95))
c50<-contourLines(gps_kde$eval.points[[1]],gps_kde$eval.points[[2]],gps_kde$estimate,level=hts[1])
c75<-contourLines(gps_kde$eval.points[[1]],gps_kde$eval.points[[2]],gps_kde$estimate,level=hts[2])
c95<-contourLines(gps_kde$eval.points[[1]],gps_kde$eval.points[[2]],gps_kde$estimate,level=hts[3])
Polyc50<-Polygon(c50[[1]][-1])
Polysc05<-Polygons(list(Polyc50),"c50")
Polyc75<-Polygon(c75[[1]][-1])
Polysc75<-Polygons(list(Polyc75),"c75")
Polyc95<-Polygon(c95[[1]][-1])
Polysc95<-Polygons(list(Polyc95),"c95")
spObj<-SpatialPolygons(list(Polysc05,Polysc75,Polysc95),1:3)
axu.df<-data.frame(ID=1:length(spObj))
rownames(axu.df)<-c("c50","c75","c95")
spDFObj<-SpatialPolygonsDataFrame(spObj, axu.df)
plot(spDFObj)

He probado muchas formas (principalmente transformando en raster y luego creando curvas de nivel) pero todas ellas han fallado porque sus resultados eran siempre significativamente diferentes al original KDE-contours-plot (plot.kde).

¿Cómo puedo exportar las curvas de nivel de mi KDE (50%, 75% y 95%) como un shapefile de polígono o polilínea para poder cargarlo en QGIS para su posterior análisis (igual que el mostrado arriba)?

Answer 1

En primer lugar, coaccione la matriz ks kde a una trama. Tenga en cuenta que tenemos que transponer la matriz kde utilizando flip .

library(raster)
library(ks)

r <- raster(extent(0,12,0,11), nrows=nrow(xykde$estimate), ncols=ncol(xykde$estimate))  
  r[] <- xykde$estimate
  r <- flip(r, direction='y')

A continuación, utilice el rasterToContour para crear un objeto de línea de contorno SpatialLinesDataFrame. Se pasa el cout vector, de nuestro objeto ks, al argumento niveles. Esto le dice a la función dónde crear los cortes del contorno.

rcont <- rasterToContour(r, levels=xykde$cont)

Aquí sólo extraemos los valores asociados a los volúmenes de contorno que nos interesan y subconjuntamos el objeto de línea sp en consecuencia. El objeto resultante contendrá contornos para los volúmenes del 50%, 75% y 95%.

( cont.values <- xykde$cont[grep(paste(c("50","75","95"), collapse = "|"), 
                            names(xykde$cont))] )

rcont.gt50 <- rcont[which(rcont@data[,1] %in% cont.values),]    

Ahora, vamos a trazar los resultados

par(mfrow=c(2,2))
  plot(xykde,display="image", main="KDE")
  plot(r, main="raster of KDE")
  plot(rcont,main="All contours (1-99%)")
  plot(rcont.gt50, main="50%, 75% and 95% volume contours")

Desde aquí se puede escribir el objeto de curvas de nivel en un shapefile utilizando rgdal::writeOGR , raster::shapefile o una de las funciones de escritura de maptools.

Answer 2

También he probado las dos respuestas anteriores con mi conjunto de datos, y puedo decir que los resultados no son los mismos. Buscando en la web, he encontrado esta solución:

Primero, convierta su objeto kde en un SpatialGridDataFrame

library(sp)

spkde <- image2Grid(list(x = xykde$eval.points[[1]], 
                         y = xykde$eval.points[[2]], 
                         z = xykde$estimate))

Asegurémonos de que tiene el mismo aspecto que un objeto kde y que el marco de datos espacial

#with kde object

image(xykde$eval.points[[1]], xykde$eval.points[[2]], xykde$estimate)

#with sp object
contour(spkde, add = TRUE)

He comprobado que la forma del kde es la misma, pero hay una ligera diferencia cuando se utiliza el "contorno" para las líneas de contorno.

Así que en este punto te recomendaría que exportaras el SpatialGridDataFrame como un raster y lo convirtieras en un polígono en QGIS directamente, pero si esto es lo suficientemente cercano para ti, puedes extraer las líneas de contorno fácilmente.

#convert into raster
r <- raster(spkde)
r.cont <- rasterToContour(r,nlevels=20)
plot(r.cont)

#Export into a shapefile
library(maptools)
writeLinesShape(r.cont,"~mycontours")

Parte de la solución vino de este enlace: https://hypatia.math.ethz.ch/pipermail/r-help/2009-November/413062.html

Answer 3

# Get polygons
hts <- contourLevels(xykde, prob = c(0.5, 0.75, 0.95))
c50 <- contourLines(xykde$eval.points[[1]], xykde$eval.points[[2]], xykde$estimate, level = hts[1])
c75 <- contourLines(xykde$eval.points[[1]], xykde$eval.points[[2]], xykde$estimate, level = hts[2])
c95 <- contourLines(xykde$eval.points[[1]], xykde$eval.points[[2]], xykde$estimate, level = hts[3])

# Convert in sp object
library(sp)
Polyc50 <- Polygon(c50[[1]][-1])
Polysc05 <- Polygons(list(Polyc50), "c50")
Polyc75 <- Polygon(c75[[1]][-1])
Polysc75 <- Polygons(list(Polyc75), "c75")
Polyc95 <- Polygon(c95[[1]][-1])
Polysc95 <- Polygons(list(Polyc95), "c95")
spObj <- SpatialPolygons(list(Polysc05,Polysc75,Polysc95), 1:3)

axu.df <- data.frame( ID=1:length(spObj))
rownames(axu.df) <- c("c50", "c75", "c95")
spDFObj <- SpatialPolygonsDataFrame(spObj, axu.df) 

# Save
library(rgdal)
writeOGR(spDFObj, dsn = '.', "levels-poly", driver="ESRI Shapefile")