4

I have shapefiles containing the location of several individuals along several months. Positions were documented with a 30 minutes interval. Tables contain the following data fields: collar (ID), date, time, lat, long, altitude.

I would like to compute the distance between each pair of points of each individual, which have been taken by 30 min or 1 hour interval apart? a number of positions and want to calculate the length between the positions - per hour, and per day.

Can someone guide the way to achieve it using R?

  • Do you mean that you want to compute the distance between each pair of points per hour / day (i.e. Speed/velocity) - or would you like to compute the velocity for each individual for each following points, sorted by time? – dof1985 Mar 25 '15 at 12:39
  • hmmm.I guess it would be for each pair of point e.g. 30 min or one hour apart in time for each individual. – Nicolai Jørgensen Mar 25 '15 at 12:41
  • What data type do you use? SHP file? .csv? etc. Also what is your data structure, e.g. variables type? – dof1985 Mar 25 '15 at 12:45
  • i use shape files with date, time, lat, long, altitude, collar (ID) – Nicolai Jørgensen Mar 25 '15 at 12:49
5

It all depends on the level of accuracy that you need.

A coarse approach would be the spherical law of cosines. This has issues with small distances - some say that it is around 1km, others say down to a few meters.

A better approach would be the Haversine Formula. This works well, however, it doesn't take into account that the earth is not actually a sphere, but a spheroid.

The "best" approach would be Vincenty's formula.

You can find a very nice implementation of all of the above at R-bloggers written by Mario Pineda-Krch.

Vincenty's formula is also implemented in R as distance() in the SDMTools-package.

A crude R-script to outline how it could be done can be found below.
It should be noted that I am by no means good at R. The script is likely inefficient and redundant. The script doesn't use the time-field in the current state, but using the lubridate-package to interpret the time and date fields will allow working with the time dimension in a decent manner.
Currently, it is assumed that the points come in the right order and that only one 'Collar ID' is present in the shapefile. These things can be dealt with by sorting the data entries by the time-field and subsetting based on the 'Collar ID'.

library(SDMTools)
library(rgeos)
library(maptools)
library(lubridate)
points <- readShapePoints("C:/Users/mlra/Desktop/FunWithDistance/DummyBird.shp")

lon <- as.vector(points@data$lon)
lat <- as.vector(points@data$lat)
z <- as.vector(points@data$z)
time <- parse_date_time(as.vector(points@data$time),"H!M!")
#The flight vector will be the distance between the current point and the previous point.
flight <- vector()
#The first element of the flight-vector will have to be zero, since no distance have been traveled
flight[1] <- 0
#We iterate from 2 to the number of observations.
for(i in 2:NROW(lon)){
  Vincenty <- distance(lat1=lat[i-1],lon1=lon[i-1],lat2=lat[i],lon2=lon[i])
  DirectDistance <- as.numeric(Vincenty[5]) # the fifth element of the output frame is the distance between the points.
  #To be more accurate, we also take into account difference in altitude between the points
  AltitudeChange <- abs(z[i]-z[i-1])
  if(AltitudeChange!=0){
    flight[i] <- sqrt(DirectDistance^2+AltitudeChange^2)
  } else {
    flight[i] <- DirectDistance
  } 
}
ResultingFrame <- data.frame(lon,lat,z,time,flight)
  • Hi. Thank you for your reply. The actual distances between the points are max 500 meters i think, but i'll look into your suggestiions :-) ... but how is it implemented in QGIS? – Nicolai Jørgensen Mar 25 '15 at 13:32
  • In that case, I'd suggest using Vincenty / the tool in SDMTools. – Mikkel Lydholm Rasmussen Mar 25 '15 at 13:34
  • and the SDMTools can be incorporated in QGIS?? – Nicolai Jørgensen Mar 25 '15 at 13:44
  • Well, strictly speaking, yes - conservationecology.wordpress.com/2013/08/14/… - but you pretty much write the script in R and then strap it onto QGIS. I'd personally stick to R for this. – Mikkel Lydholm Rasmussen Mar 25 '15 at 13:50
  • Ok. Thanks. Then i've got two things to do. Writing the script and find out how to strap it into QGIS... – Nicolai Jørgensen Mar 25 '15 at 13:52
3

After a little editing i used this script and its working fine. Now i just have to do it for each ID and for each day :-)

library(SDMTools)
library(rgeos)
library(maptools)
library(lubridate)

#points <- readShapePoints("C:/Users/mlra/Desktop/FunWithDistance/DummyBird.shp")
points <- tellus_2012_2013

lon <- as.vector(points$long)
lat <- as.vector(points$lat)
z <- as.vector(points$alt)
time <- parse_date_time(as.vector(points$time),"%H:%M:%S")

#The flight vector will be the distance between the current point and the previous point.
flight <- vector()

#The first element of the flight-vector will have to be zero, since no distance have been traveled
flight[1] <- 0

#We iterate from 2 to the number of observations.
for(i in 2:NROW(lon)){
  Vincenty <- distance(lat1=lat[i-1],lon1=lon[i-1],lat2=lat[i],lon2=lon[i])
  DirectDistance <- as.numeric(Vincenty[5]) # the fifth element of the output frame is the distance between the points.
  #To be more accurate, we also take into account difference in altitude between the points
  AltitudeChange <- abs(z[i]-z[i-1])
  if(AltitudeChange!=0){
    flight[i] <- sqrt(DirectDistance^2+AltitudeChange^2)
  } else {
    flight[i] <- DirectDistance
  } 
}

 ResultingFrame <- data.frame(lon,lat,z,time,flight)
  • I can get ID and date appended to the resulting frame and might be able to do something in a e.g. excel to do some sorting... – Nicolai Jørgensen Mar 26 '15 at 10:40
  • Would it be possible to make a "for each loop" for each id, for each year, for each day, for each id and for each something else in the above script? – Nicolai Jørgensen Mar 27 '15 at 6:44
  • Anyone for the QGIS bit?? – Nicolai Jørgensen Mar 27 '15 at 9:49
  • For the QGIS bit you should ask another question that focuses on that. – PolyGeo Oct 29 '17 at 6:32

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.