Kernel Density Estimation Bootstrap

[Christian Gredzens]

Hi All,

I would like to run a bootstrap within R on a spatial dataset using a kernel density estimation function. Can anyone help me with this? I am new to R and thus have very little experience. It seems relatively easy to run the actual bootstrap but I am having problems in the following areas:

1) I don't know how to import spatial data into R and assign it an appropriate projection
2) I need to setup a kernel function using a specific bandwidth (likelihood cross-validation (CVh)) and create the script so it will run within the bootstrap script
3) I need to export the bootstrap data as an ArcGIS 10 file.

Best Regards,
Christian Gredzens

Christian Gredzens
MSc Student
School of Marine and Tropical Biology
James Cook University
Townsville, QLD, Australia 4811
Christian...@my.jcu.edu.au
Phone: 04 3875 3652

[Stewart Macdonald]

Hi Christian,

What sort of data are you trying to import? A series of lat/longs that are in a text (or Excel) file? If so, the rGDAL library can be used to create spatial points:

###############
library(rgdal)
# if you don't already have this installed, you can do: install.packages('rgdal')

# read in the data from a tab-delimited text file that has a header line
rawPoints <- read.table('/data.csv', header=TRUE, sep='\t')

# format as a SpatialPoint class and define the projection/datum (WGS84)
q <- cbind(rawPoints$recordLon, rawPoints$recordLat)
colnames(q) <- c("x", "y")
qqq <- SpatialPoints(q, proj4string=CRS("+init=epsg:4326"))
###############

Note: the above is a copied-and-pasted fragment from a working script of mine. I assume this fragment will work for you.

I can't help you with your second question, and the answer to your third question will depend on the format of the bootstrap data (e.g., more points, polygons, raster).


Hope this helps,

Stewart

> --
> An R group for questions, tips and tricks relevant to spatial ecology and climate change.
> All R questions welcome.
> ---
> You received this message because you are subscribed to the Google Groups "Tropical R" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to tropical-r+...@googlegroups.com.
> To post to this group, send an email to tropi...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msg/tropical-r/-/DgXQECjrCaoJ.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Leila Brook]

Hi Christian,
There is a package called adehabitatHR that does kernel density estimation with LSCV bandwidth using a function called kernelUD(). I don't know if this will is similar to what you want? You will need to have the data in UTM and in SpatialPoints or SpatialPointsDataFrame format, which you can do using the sp package.

There is a great vignette for adehabitatHR that will help with data formatting and running the KDE. I exported the utilisation distribution for ArcGIS using writeRaster() in the raster package, but you can also export as contours.

Cheers
Leila

________________________________________
From: tropi...@googlegroups.com [tropi...@googlegroups.com] on behalf of Stewart Macdonald [stewartm...@gmail.com]
Sent: Friday, 3 May 2013 5:12 PM
To: tropi...@googlegroups.com
Subject: Re: [Trop-R] Kernel Density Estimation Bootstrap

[Christian Gredzens]

Hi Stewart,

Thanks for the help. I'm using a .csv file with lats/longs in UTM entered as:

x	y
648501	8910400
641686	8917264
645984	8914328
647058	8913794
649659	8912925
649613	8912828

....            .....

I can get R to read the file using rgdal and get the following:

> test<-read.table("test.csv")
> test
                V1
1              x,y
2   648501,8910400
3   641686,8917264
4   645984,8914328
5   647058,8913794
6   649659,8912925
7   649613,8912828
.....

But when I do the next step I get the following output:

> q<-cbind(test$recordLon, test$recordLat)
> q
NULL
> colnames(q)<-c("x","y")
Error in `colnames<-`(`*tmp*`, value = c("x", "y")) :
  attempt to set colnames on object with less than two dimensions


I'm sure this is just my inexperience with R, but what am I missing?

I tried it this way as well with the same result:

> q<-cbind(test$x, test$y)
> q
       [,1]
  [1,]  268
  [2,]  111
  [3,]  221
...
> colnames(q)<-c("x","y")
Error in `colnames<-`(`*tmp*`, value = c("x", "y")) :
  length of 'dimnames' [2] not equal to array extent


Also, once I run the bootstrap I would like to export the results as a polygon layer.

Cheers,
Christian

[Christian Gredzens]

Leila,

Thanks for your reply. I'm using likelihood cross validation (CVh) instead of least squares cross validation (LSCV). To my knowledge adehabitatHR doesn't run CVh so I think I will need to write a script or find another program which can calculate CVh.

I used GME from spatialecology.com to calculate all of my KDEs, but I don't think GME runs bootstrapping. So I'm trying to run it through the bootstrapping package in R. Here is what I need to do:

1) Create my spatial points layer
=mydata

2) Run the bootstrapping script:

>mydataboot<-boot(mydata, (KDE script), R=1000)

mydata= my spatial points file
(KDE script)= what formula to run the bootstrap on
R=1000 = the number of iterations I would like to run

The issue I'm having is getting the bootstrapping script to run the KDE formula=(KDE script)

Cheers,
Christian

[Stewart Macdonald]

Looks like you need to specify the field delimiter:

> test<-read.table("test.csv", sep=",")

This tells R that your fields are separated by commas. At the moment it's just treating the entire line as the sole field.

Also, if your column names are already x and y, you can remove the 'colnames(q)<-c("x","y")' line, and alter the other line thusly:

q <- cbind(test$x, test$y)



Stewart

> To view this discussion on the web, visit https://groups.google.com/d/msg/tropical-r/-/8bjSxroW9XMJ.

[Lauren Hodgson]

Hi Christian,

This formatting shows that the data has been read in as a single column with 'V1' as the column name and 'x,y' as the value of row 1, and so on:

> test<-read.table("test.csv")
> test
                V1
1              x,y
2   648501,8910400
3   641686,8917264
4   645984,8914328
5   647058,8913794
6   649659,8912925
7   649613,8912828
.....

You can check how many columns there are in a data frame using 'ncol':

ncol(test)

To troubleshoot, I would first try using read.csv:

test=read.csv(test,as.is=TRUE)

ncol(test)

If this does not give you two columns, then I would check that the file really is in standard .csv format.

Lauren

To view this discussion on the web, visit https://groups.google.com/d/msg/tropical-r/-/T2kO_NrGe3MJ.

[Alice C. Hughes]

Hi All

This should be simple-but ArcGis keeps getting stuck at 22%...

I have a large shapefile, which contains multiple overlapping polygons for various species-and I would like to count for any area how many polygons overlap.

A script would be really fantastic!

Thanks in advance

Alice

Alice C. Hughes

Associate Professor

Centre for Integrative Conservation,

Xishuangbanna Tropical Botanical Garden,

Chinese Academy of Sciences

Menglun, Mengla, Yunnan 666303, P.R. China

Ph: 15198676559

achu...@xtbg.ac.cn

---

Date: Mon, 6 May 2013 14:33:29 +1000

Subject: Re: [Trop-R] Kernel Density Estimation Bootstrap

From: lhodg...@gmail.com
To: tropi...@googlegroups.com

[Alice C. Hughes]

Hi All

Could anyone help me with a script which would scroll through a list of polygons of species ranges, then go through either an excel, or point file and list how many points for each species fall inside, and how many fall outside that species polygon in a table?

Thanks in advance

Alice

Alice C. Hughes

Associate Professor

Centre for Integrative Conservation,

Xishuangbanna Tropical Botanical Garden,

Chinese Academy of Sciences

Menglun, Mengla, Yunnan 666303, P.R. China

Ph: 15198676559

achu...@xtbg.ac.cn

---

Date: Mon, 6 May 2013 14:33:29 +1000

Subject: Re: [Trop-R] Kernel Density Estimation Bootstrap

From: lhodg...@gmail.com
To: tropi...@googlegroups.com