4

I'm performing a kernel (KDE) home range analysis of an individual with the package 'ks' in R. I'd like to avoid using Geospatial Modelling Environment (GME) for implementing the output in ArcGIS, because I'm using QGIS which is not compatible with GME. I'd like run the analysis in R (see example below) and import the output as a shapefile into 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))

original KDE-contours

These contours are exactly what I want to export as a shapfile, but there is still no reasonable possibility described in the internet to do this.

Adapting @Jeffrey Evans' Formula resulted in a KDE-plot which differed significantly from the original plot.kde function.

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")

@Jeffrey Evans' plots

Further analyzes are not possible with this result. Also the method of @Juan Antonio Roldán Díaz resulted in non-applicable results.

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)

@Juan Antonio Roldán Díaz' plot

I have tried many ways (mainly tranforming into raster and then creating contour lines) but they all failed because their results were always significantly different to the original KDE-contours-plot (plot.kde).

How do I export the contour lines of my KDE (50%, 75% and 95%) as a polygon or polyline shapefile so I can load it into QGIS for further analysis (same as shown above)?

2

I have also tried both answers above with my dataset, and I can relate that results are not the same. By searching on the web, I found this solution:

First, convert your kde object into a SpatialGridDataFrame

library(sp)

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

Let's make sure it looks exactly the same as both a kde object and the spatial dataframe

#with kde object

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

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

I have found that the shape of the kde is the same, but there is a slight difference when you use the 'contour' for the contour lines.

So at this point I would recommand that you export the SpatialGridDataFrame as a raster and convert into a polygon on QGIS directly, but if this is close enough for you, you can extract the contour lines easily.

#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")

Part of the solution came from this link: https://hypatia.math.ethz.ch/pipermail/r-help/2009-November/413062.html

2

First coerce the ks kde matrix to a raster. Please note that we have to transpose the kde matrix using 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')

Then, use the rasterToContour to create a SpatialLinesDataFrame contour line object. You pass the cout vector, from our ks object, to the levels argument. This tells the function where to create contour breaks.

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

Here we just extract the values associated with the contour volumes that you are interested in and subset the sp line object accordingly. The resulting object will contain contours for the 50%, 75%, and 95% volumes.

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

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

Now, let's plot the results

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")

From here you can write the contour line object to a shapefile using rgdal::writeOGR, raster::shapefile or one of the write functions in maptools.

4
  • You're missing library(raster) anywhere! – Spacedman May 25 '18 at 17:13
  • 1
    Does rasterToContour(r, levels=cont.values) get the required contours without having to compute all the contours and then grepping out the required ones? – Spacedman May 25 '18 at 17:17
  • Yes, not sure why I went the long way. The most efficient way would be: rasterToContour(r, levels=contourLevels(xykde, prob = c(0.5, 0.75, 0.95))) – Jeffrey Evans May 25 '18 at 18:00
  • Thanks for your response! Unfortunately, when perfoming KDE on real GPS-data, the results of your formula differ significantly from the results of the plot.kde function. More informations in my edited question. – Raoul Reding Sep 1 '19 at 18:13
0
# 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")
3
  • maptools::ContourLines2SLDF will do most of that for you. – Spacedman May 25 '18 at 17:10
  • ...and you can't rely on contour lines being polygons - they can stop at the edges of the region. – Spacedman May 25 '18 at 17:11
  • Thanks for your response! Unfortunately, when perfoming KDE on real GPS-data, the results of your formula differ significantly from the results of the plot.kde function. More informations in my edited question. – Raoul Reding Sep 1 '19 at 18:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.