[R-sig-eco] adehabitat and Lat Long data
Pinaud David
pinaud at cebc.cnrs.fr
Tue Aug 28 09:37:13 CEST 2012
Hi Mike,
you could use these homemade First Passage Time functions, based on
Fauchald & Tveraa, Ecology (2003), 84:282-288 and the trip package. It
uses distance calculation in km if longlat=T (DD) or in CRS unit if data
are projected. You could use then Lavielle's method, used by Barraquand
and Benhamou (Ecology, 2008).
##### an example:
source("FPT4.r")
library(trip)
library(adehabitatLT)
### FPT analysis:
## use the FPT4.r function collection. Data are a trip object ("tr2")
scales <- c(seq(0.01, 0.1, by=0.01), seq(0.2, 1, by=0.1), seq(1.2, 5,
by=0.2), 6:30) # the scales at which FPT (=radius) is calculated
tr2.cut <- FPT.cut(tr2, d=0.05) # cut the trip at every d distance (=
add intermediate locs)
loc.fpt <- FPT.calc(tr2.cut, r=scales) # calculate FPT, can be quite
long (several minutes)...
# draw the log(FPT)=f(scales) plot
savar <- FPT.vargraph(loc.fpt, r=scales, name="tr2", plotit=T, saveit=F)
# peaks of variance at 0.1 and 8 km -> ARS at these scales
# draw FPT at ARS scales (here 0.1 km) against
time since departure:
plot(x=getTimeID(loc.fpt)[, 1], y=loc.fpt at data[, "r0.1"], t="l",
main="FPT at 0.1km scale", xlab="Time since departure", ylab="FPT at 0.1
km")
# to find where (when) ARS occured. This segmentation is based
on Lavielle's method, used by Barraquand and Benhamou (Ecology, 2008).
# ex at 0.1 km
l.1 <- lavielle(loc.fpt at data[!is.na(loc.fpt at data$r0.1), "r0.1"], Lmin=2,
Kmax=20, ld=2)
chooseseg(l.1) # give the number of segments to keep, here we should use
K=3,
fp1 <- findpath(l.1, 3) # cut the path into K segments and give the
row ID of segement limits
HTH
David
Le 27/08/2012 16:37, Mike a écrit :
> Hi
>
> I have ARGOS satellite tag data from some some very long distance traveling
> sharks. These
> critters cover large chunks of the Atlantic; coordinates are in
> decimal degrees.
>
> It would love to apply the residenceTime function in the adehabitat package
> to try and segment the paths into broad homogeneous bouts of movement (i.e.
> migration, seasonal residence/feeding areas etc). But as far as I can tell
> from the documentation and the examples it seems that the data must be in
> some sort of Cartesian coordinate system. Conversion to UTM's wouldn't be
> practical in
> this case because the movement tracks encompass many UTM zones. Is there
> a way to do the residence time analysis in adehabitat using data in decimal
> degrees that perhaps I am just missing?
>
> If not, does anybody have any code or recommendations for another way to go
> about this? Any insight appreciated!
>
> Thanks
> -Mike
>
>
>
>
> --
> View this message in context: http://r-sig-ecology.471788.n2.nabble.com/adehabitat-and-Lat-Long-data-tp7577571.html
> Sent from the r-sig-ecology mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
--
***************************************************
Dr. David PINAUD
Ingénieur de Recherche "Analyses spatiales"
Centre d'Etudes Biologiques de Chizé - CNRS UPR1934
79360 Villiers-en-Bois, France
poste 485
Tel: +33 (0)5.49.09.35.58
Fax: +33 (0)5.49.09.65.26
http://www.cebc.cnrs.fr/
***************************************************
-------------- next part --------------
########### FIRST-PASSAGE TIME ANALYSIS ###########
# v.25/06/2012
# David PINAUD, pinaud at cebc.cnrs.fr, June 2012. Developped with R 2.14.2 (2012-02-29)
# Many thanks to LôRan Dubroca and Pascal Monestier
# Based on Fauchald & Tveraa, Ecology (2003), 84:282-288
#
# - require "sp", "fields" & "trip" libraries
# - input file: a trip object
#
# Analysis in 3 steps:
# - FPT.cut() calculates intermediate fixes (X, Y and date) for each successive distance d
# - FPT.calc() calculates First-Passage Time for each intermediate fix and for each spatial scale "r"
# - FPT.vargraph() calculates the var(log(fpt)) and plot it for each spatial scale r.
#
#
## FPT.cut calculates intermediate fixes (X, Y and date) for each successive distance d (in km if longlat=T or in CRS unit), from a trip object.
## Returns a file to be passed to FPT.calc() analysis. Possibility to keep "extra data".
## Version vectorielle + rapide.....
FPT.cut <- function(tr, d, longlat=!is.projected(tr))
{
require(trip)
tor <- getTORnames(tr) # get the DateTime and ID
coordnames <- dimnames(coordinates(tr))[[2]] # get the coord names
dist01 <- c(0, trackDistance(coordinates(tr), longlat = longlat)) # distance between subsequent locations, in km if longlat=T or in d (CRS) unit
Xv <- c(0, coordinates(tr)[2:nrow(coordinates(tr)), 1] - coordinates(tr)[1:(nrow(coordinates(tr))-1), 1]) # longitude vectors
Yv <- c(0, coordinates(tr)[2:nrow(coordinates(tr)), 2] - coordinates(tr)[1:(nrow(coordinates(tr))-1), 2]) # latitude vectors
lag <- dist01 / d # number of intermediate fixes according to d /!\ d should be in CRS unit or in km if longlat=T
rm(dist01)
nbfix <- trunc(lag)
nbfix <- ifelse(((nbfix == lag) & (lag > 0)), lag-1, nbfix ) # case where dist is a multiple of d and lag==0
pos0 <- cumsum(nbfix+1) # index of locations
mc <- as.data.frame(matrix(data=FALSE, nr=(sum(nbfix)+length(nbfix)), nc=5))
names(mc) <- c(tor, coordnames, "Origin")
mc[pos0, 1:4] <- data.frame(getTimeID(tr), coordinates(tr))
mc[pos0, 5] <- TRUE
Xvr <- rep(Xv/lag, nbfix+1) # vectors X
Yvr <- rep(Yv/lag, nbfix+1) # vectors Y
rm(Yv, Xv)
Dvr <- c(0, rep((getTimeID(tr)[2:nrow(getTimeID(tr)), 1] - getTimeID(tr)[1:(nrow(getTimeID(tr)) - 1), 1]) / lag[2:nrow(getTimeID(tr))], (nbfix+1)[2:length(nbfix)])) # vectors Date
seqr <- sequence(nbfix+1) # sequence of multipliers
Xb <- c(rep(coordinates(tr)[1:(length(nbfix)-1), 1], (nbfix+1)[2:length(nbfix)]), 0) # basis X
Yb <- c(rep(coordinates(tr)[1:(length(nbfix)-1), 2], (nbfix+1)[2:length(nbfix)]), 0) # basis Y
Dtb <- c(rep(getTimeID(tr)[1:(length(nbfix)-1), 1], (nbfix+1)[2:length(nbfix)]), 0) # basis DateTime
rm(lag)
# calculus:
mc[mc$Origin==FALSE, coordnames[1]] <- Xb[mc$Origin==FALSE] + Xvr[mc$Origin==FALSE] * seqr[mc$Origin==FALSE]
mc[mc$Origin==FALSE, coordnames[2]] <- Yb[mc$Origin==FALSE] + Yvr[mc$Origin==FALSE] * seqr[mc$Origin==FALSE]
mc[mc$Origin==FALSE, tor[1]] <- Dtb[mc$Origin==FALSE] + Dvr[mc$Origin==FALSE] * seqr[mc$Origin==FALSE]
rm(Xvr, Yvr, Dvr, seqr, Xb, Yb, Dtb)
mc[, tor[1]] <- mc[, tor[1]] + as.POSIXct("1970-01-01 00:00:00", "GMT")
mc[, tor[1]] <- as.POSIXct(mc[, tor[1]], tz="GMT")
mc[, tor[2]] <- getTimeID(tr)[1, 2]
# extra data
if(ncol(tr at data)>4) {
"%w/o%" <- function(x, y) x[!x %in% y] #-- x without y
cd <- tr at data[, colnames(tr at data) %w/o% tor]
kept <- as.data.frame(matrix(data=NA, nr=(sum(nbfix)+length(nbfix)), nc=ncol(cd)))
names(kept) <- colnames(cd)
kept[pos0, ] <- cd
mc <- data.frame(mc, kept)
}
# trip
coordinates(mc) <- coordnames # is now a 'sp' object (SpatialPointsDataFrame)
proj4string(mc) <- CRS(proj4string(tr)) # set the coord system / projection
mc <- trip(mc, tor) # a trip object, "locs$device_info_serial" as factor
return(mc)
rm(nbfix, pos0, mc, kept, cd)
}
## "FPT.calc" calculates First-Passage Time for each successive fix of "tr" object (coming from FPT.cut() ) and for each spatial scale "radii".
## "tr" is the result of "FPT.cut", "radii" is a sequence of each spatial scale to be tested in km (if tr is longlat) or in coordinate unit si tr is projected.
FPT.calc <- function (tr, radii)
{
options(digits = 10)
startTime <- Sys.time()
fpt2 <- NULL
n <- nrow(tr at data)
cat("0%")
for (i in 1:n)
{
disttemp <- spDistsN1(pts=as(tr, "SpatialPointsDataFrame")[i:n,], pt=as(tr, "SpatialPointsDataFrame")[i,], longlat = !is.projected(tr))
distclas <- cut(disttemp, breaks=c(0, radii, (max(radii)+ceiling(diff(range(radii))/length(radii)))), include.lowest = T)
fpt.fwd <- tapply(getTimeID(tr)[i:n, 1], distclas, min) - as.numeric(getTimeID(tr)[i, 1])
disttemp <- spDistsN1(pts=as(tr, "SpatialPointsDataFrame")[1:i,], pt=as(tr, "SpatialPointsDataFrame")[i,], longlat = !is.projected(tr))
distclas <- cut(disttemp, breaks=c(0, radii, (max(radii)+ceiling(diff(range(radii))/length(radii)))), include.lowest = T)
fpt.bwd <- as.numeric(getTimeID(tr)[i, 1]) - tapply(getTimeID(tr)[1:i, 1], distclas, max)
fpttotal <- apply(rbind(fpt.bwd, fpt.fwd), 2, sum)
fpt2 <- rbind(fpt2, fpttotal)
# time bar:
tb <- cbind(trunc(n/10)*1:10, 1:10*10)
if(any(i==tb[,1])) cat(paste("|| ", tb[i==tb[,1],2], "%, ", round(c(difftime(Sys.time(), startTime, unit="mins"))), "min", sep=""))
}
fpt2 <- as.data.frame(fpt2)
names(fpt2) <- paste("r", c(0, radii), sep="")
fpt2$r0 <- NULL
tr at data <- cbind(tr at data, fpt2)
return(tr)
}
## FPT.vargraph calculates the var(log(fpt)) for a file from FPT.calc()and plot it for each radius.
## Return a data.frame and a plot for the individual "name". Plot can be saved.
FPT.vargraph <- function(tr.fpt, radii, name="", modes.method=1, plotit=F, saveit=F, type=c("wmf", "png", "jpeg", "jpg", "bmp", "ps", "pdf"))
{
r.ind <- radii
r.char <- paste("r",radii, sep="")
r.ind <- as.data.frame(t(radii))
names(r.ind) <- r.char
r.ind[1,] <- NA
fptr <- subset(tr.fpt at data, select=r.char)
for (j in 1:(length(radii)))
{
jr<-fptr[, j]
if(length(na.omit(log(jr))) > 1) { r.ind[1, j] <- stats::var(log(jr), na.rm=T) }
}
rfpt <- radii[r.ind == max(r.ind, na.rm=T)]
if (plotit)
{
if(modes.method==2)
{
pm <- posmodes2(radii, r.ind, numbermodes2(radii, r.ind))
} else {
pm <- posmodes(radii, r.ind, numbermodes(radii, r.ind))
}
plot(radii, r.ind, t="l", xlab="Spatial scale r (km)", ylab="S(log(fpt))", xlim=c(0, max(radii)),
ylim=c(0, max(r.ind, na.rm=T)), main=paste(name, ", FPT analysis", sep=""))
abline(v=radii, col="grey", lty=3)
points(radii, r.ind, t="l")
abline(v=rfpt, col="red", lty=2)
if(length(pm)>0) text(x=pm, y=rep(c(0, max(r.ind, na.rm=T)/20), length.out=length(pm)), labels=as.character(pm), cex=0.85, font=4)
}
if(plotit & saveit)
{
savePlot(filename=paste(name,"vargraf", sep=""), type=type)
}
return(r.ind)
}
### FPT.cutsim calculates intermediate fixes (X, Y and date) for each successive distance d (in km) for simulated data (cartesian coordinates).
### Returns a file to be passed to FPT.calcsim() analysis.
#FPT.cutsim <- function (file, Long, Lat, Date, d, datakept=NULL)
#{
# options(digits=10)
# require(fields)
# a <- data.frame(Date, Long, Lat, Origin=rep("FALSE", length(Long)))
#
# coord <- matrix(c(a$Long, a$Lat), length(a$Long), 2) # XY matrix, to calculate distance
# coord <- coord[-length(coord[, 1]), ]
# coord2 <- matrix(c(a$Long[2:length(a$Long)], a$Lat[2:length(a$Lat)]), length(a$Long)-1, 2)
# dist01 <- c(0, sqrt((coord[,1]-coord2[,1])^2 + (coord[,2]-coord2[,2])^2)) # dist to previous fix
# rm(coord, coord2)
#
# Xv <- c(0, a$Long[2:length(a$Long)] - a$Long[1:(length(a$Long) - 1)]) # longitude vectors
# Yv <- c(0, a$Lat[2:length(a$Lat)] - a$Lat[1:(length(a$Lat) - 1)]) # latitude vectors
#
# lag <- dist01 / d # number of intermediate fixes according to d
# rm(dist01)
# nbfix <- trunc(lag)
# nbfix <- ifelse (((nbfix == lag) & (lag > 0)), lag-1, nbfix ) # case where dist is a multiple of d and lag==0
#
# pos0 <- cumsum(nbfix+1) # position des points GPS
#
# mc <- as.data.frame(matrix(data=FALSE, nr=(sum(nbfix)+length(nbfix)), nc=4))
# names(mc) <- c("Date", "Long", "Lat", "Origin")
# mc[pos0, 1:3] <- a[, 1:3]
# mc[pos0, 4] <- TRUE
#
# Xvr <- rep(Xv/lag, nbfix+1) # vecteurs X
# Yvr <- rep(Yv/lag, nbfix+1) # vecteurs Y
# rm(Yv, Xv)
# Dvr <- c(0, rep((a$Date[2:length(a[,1])]-a$Date[1:(length(a[,1])-1)])/lag[2:length(a[,1])],(nbfix+1)[2:length(nbfix)])) # vecteurs Date
# seqr <- sequence(nbfix+1) # sequence des multiplicateurs
# Xb <- c(rep(a$Long[1:(length(nbfix)-1)], (nbfix+1)[2:length(nbfix)]), 0) # bases X
# Yb <- c(rep(a$Lat[1:(length(nbfix)-1)], (nbfix+1)[2:length(nbfix)]), 0) # bases Y
# Db <- c(rep(a$Date[1:(length(nbfix)-1)], (nbfix+1)[2:length(nbfix)]), 0) # bases Date
# rm(a, lag)
#
# # calcul:
# mc$Date[mc$Origin==FALSE] <- Db[mc$Origin==FALSE] + Dvr[mc$Origin==FALSE]*seqr[mc$Origin==FALSE]
# mc$Long[mc$Origin==FALSE] <- Xb[mc$Origin==FALSE] + Xvr[mc$Origin==FALSE]*seqr[mc$Origin==FALSE]
# mc$Lat[mc$Origin==FALSE] <- Yb[mc$Origin==FALSE] + Yvr[mc$Origin==FALSE]*seqr[mc$Origin==FALSE]
# rm(Xvr, Yvr, Dvr, seqr, Xb, Yb, Db)
#
# # ajout de datakept
# if(!is.null(datakept)) {
# kept <- as.data.frame(matrix(data=NA, nr=(sum(nbfix)+length(nbfix)), nc=ncol(datakept)))
# names(kept) <- names(datakept)
# kept[pos0, ] <- datakept
# mc <- data.frame(mc, kept)
# }
#
# return(mc)
# rm(nbfix, pos0, mc, kept)
# }
#
### "FPT.calcsim" calculates First-Passage Time for each successive fix of "a" and for each spatial scale "r" for simulated data.
## "a" is the result of "FPT.cutsim", "r" is a sequence of each spatial scale to be tested in km.
#FPT.calcsim <- function (a, r)
#{
# options(digits = 13)
# require(fields)
# fpt2 <- NULL
# n <- nrow(a)
# for (i in 1:n)
# {
# disttemp <- rdist(matrix(c(a$Lat[i], a$Long[i]), length(a$Long[i]), 2), matrix(c(a$Lat[i:n], a$Long[i:n]), length(a$Long[i:n]), 2))
# distclas <- cut(disttemp, breaks=c(0, r, (max(r)+ceiling(diff(range(r))/length(r)))), include.lowest = T)
# fptav <- tapply(a$Date[i:n], distclas, min) - a$Date[i]
# disttemp <- rdist(matrix(c(a$Lat[i], a$Long[i]), length(a$Long[i]), 2), matrix(c(a$Lat[1:i], a$Long[1:i]), length(a$Long[1:i]), 2))
# distclas <- cut(disttemp, breaks=c(0, r, (max(r)+ceiling(diff(range(r))/length(r)))), include.lowest = T)
# fptar <- a$Date[i] - tapply(a$Date[1:i], distclas, max)
# fpttotal <- apply(rbind(fptar, fptav), 2, sum)
# fpt2 <- rbind(fpt2, fpttotal)
#
# if(i/100-trunc(i/100)==0) print(paste(i,"th fix on", n))
# }
# row.names(fpt2) <- 1:n
# fpt2 <- as.data.frame(fpt2[, -1])
# names(fpt2) <- paste("r",r, sep="")
# a <- cbind(a, fpt2)
# return(a)
# rm(a)
#}
#
# "numbermodes" finds number of modes- Bradley Efron's method
# modified from code (function nmodes) supplied by Rob Tibshirini
# to return only number of modes.
numbermodes <- function(x, y)
{
n <- length(x)
mcount <- 0
for(i in 1:(n-2))
{
if(is.na(y[i +2])) break
if((y[i + 1] > y[i]) & (y[i + 1] > y[i +2]))
{
mcount <- mcount + 1
}
}
return(mcount)
}
numbermodes2 <- function(x, y) # fenêtre plus large (sur -2 et +2 points)
{
n <- length(x)
mcount <- 0
for(i in 2:(n-3))
{
if((y[i + 1] > y[i]) & (y[i + 1] > y[i +2]) & (y[i ] > y[i-1]) & (y[i + 2] > y[i +3]))
{
mcount <- mcount + 1
}
}
return(mcount)
}
# "posmodes" renvoie la position (x) des modes de "y"
# selon le nombre donnné par "numbermodes".
posmodes <- function(x, y, nmodes)
{
n <- length(x)
mcount <- 0
pos <- vector(mode="numeric", length=nmodes)
for(i in 1:(n-2))
{
if(is.na(y[i +2])) break
if((y[i + 1] > y[i]) & (y[i + 1] > y[i +2]))
{
mcount <- mcount + 1
pos[mcount] <- x[i+1]
}
}
return(pos)
}
# "posmodes2" renvoie la position (x) des modes de "y"
# selon le nombre donnné par "numbermodes2", fenêtre plus large (sur -2 et +2 points).
posmodes2 <- function(x, y, nmodes)
{
n <- length(x)
mcount <- 0
pos <- vector(mode="numeric", length=nmodes)
for(i in 2:(n-3))
{
if(is.na(y[i + 3])) break
if((y[i + 1] > y[i]) & (y[i + 1] > y[i +2]) & (y[i ] > y[i-1]) & (y[i + 2] > y[i +3]))
{
mcount <- mcount + 1
pos[mcount] <- x[i+1]
}
}
return(pos)
}
### rdist.earth2 calculates distance between loc1[i] and loc2[i].
# Rlat adjusts earth's radius to the latitude of the first fix of loc1.
rdist.earth2 <- function (loc1, loc2, Rlat = TRUE)
{
R <- ifelse(Rlat, 6378137 * (1 - 0.0033493*sin(loc1[1,2])^2)/1000, 6378.388)
if (missing(loc2))
loc2 <- loc1
coslat1 <- cos((loc1[, 2] * pi)/180)
sinlat1 <- sin((loc1[, 2] * pi)/180)
coslon1 <- cos((loc1[, 1] * pi)/180)
sinlon1 <- sin((loc1[, 1] * pi)/180)
coslat2 <- cos((loc2[, 2] * pi)/180)
sinlat2 <- sin((loc2[, 2] * pi)/180)
coslon2 <- cos((loc2[, 1] * pi)/180)
sinlon2 <- sin((loc2[, 1] * pi)/180)
pp <- (coslat1 * coslon1 * coslat2 * coslon2) + (coslat1 * sinlon1 * coslat2 * sinlon2) + (sinlat1 * sinlat2)
R * acos(ifelse(pp > 1, 1, pp))
}
# calcul de distance sur le globe terrestre (je fais confiance!!)
disttraj <- function(LAT1, LAT2, LONG1, LONG2)
{
dist1 <- (((((((LAT2-LAT1)*3600)*6335508)*((1-0.00672267*((sin(((LAT2+LAT1)/2)*2*pi/360))^2))^(-2/3))*(1/206265))^2)+(((((LONG2-LONG1)*3600)*6378388)*((1-0.00672267*(sin(((LAT2+LAT1)/2)*2*pi/360)))^(-1/2))*(1/206265)*(cos(((LAT2+LAT1)/2)*2*pi/360)))^2))^(1/2))/1000
return(dist1)
}
# calcul du ratio Lon/Lat pour les graphs (à passer dans asp=)
asp.rect <- function(Long, Lat) {
coord.m <- data.frame(mean(Long), mean(Lat))
d1.long <- rdist.earth2(coord.m, coord.m+c(1,0))
d1.lat <- rdist.earth2(coord.m, coord.m+c(0,1))
rap <- d1.lat/d1.long
return(rap)
}
# fonction distp en version longitude latitude
distp2 <- function(latp, longp, lat, long)
{
# latp et longp sont de dimension1, lat et long sont de meme dimension n
n <- length (lat)
dist <- rep(0, n)
dist <- disttraj(rep(latp, n), lat, rep(longp, n), long)
return(dist)
}
More information about the R-sig-ecology
mailing list