20 Nov

Get angle/orientation/direction of GPS track using R

Hey,

Today a bit of trSIGonometry with this function to compute the angle/orientation/direction between 2 points of a GPX file.

############################
#Function
############################
 
angle <- function(pt1, pt2) {
  conv = 180 / pi # conversion radians / degrees
  diff <- pt2 - pt1
  diff.abs <- abs(diff)
  if (diff[1]>=0 & diff[2]>=0){angle <- pi/2 - atan(diff.abs[2] / diff.abs[1])}
  if (diff[1]>=0 & diff[2]<=0){angle <- pi/2 + atan(diff.abs[2] / diff.abs[1])}
  if (diff[1]<=0 & diff[2]>=0){angle <- 3*pi/2 + atan(diff.abs[2] / diff.abs[1])}
  if (diff[1]<=0 & diff[2]<=0){angle <- pi + (pi/2- atan(diff.abs[2] / diff.abs[1]))}
  return(angle * conv)
}
 
############################
#Script
############################
 
library(rgdal)
library(raster)
library(data.table)
 
#Set folder
setwd("/home/GuruBlard/gpx/")
 
# Load file names
filenames <- list.files(pattern=".gpx$")
 
# Iteration through all gpx files
for (filename in filenames){
  # Load the GPX
 
  #ogrListLayers("export.gpx")
  gpx.trackpoints <- readOGR(dsn = filename,  layer= "track_points", stringsAsFactors = F)
 
  # Extract the coordinates as a matrix:
  gpx.tp <- coordinates(gpx.trackpoints)
 
  # #Loop through and add angle - the difference between the of the current point and the following one:
  for (tp in 1:(nrow(gpx.tp)-1)) {
    gpx.trackpoints$angle[tp] <- angle(gpx.tp[tp,], gpx.tp[tp+1,])
  }
 
  gpx.trackpoints@data$Date <- substr(as.character(gpx.trackpoints@data$time), 1,10)
  gpx.trackpoints@data$Hour <- substr(as.character(gpx.trackpoints@data$time),12,19)
  gpx.trackpoints@data$filename <- filename
  writeOGR(gpx.trackpoints, paste0(tools::file_path_sans_ext(filename), ".shp"), layer=tools::file_path_sans_ext(filename), driver="ESRI Shapefile")
}

Screenshot from 2017-11-20 14-53-09

13 Jan

Plot points with a google map background

Good Morning Followers!

Yesterday, I found a nice post here by Joseph Rickert. I produces nice plots with points and a google map background. Try it and have fun!

library(ggmap)
library(plyr)
library(gridExtra)
temp <- tempfile()
download.file("http://www.plantsciences.ucdavis.edu/plant/data.zip", temp)
connection <- unz(temp, "Data/Set3/Set3data.csv")
rice <- read.csv(connection)
names(rice) <- tolower(names(rice))
# Create a custom soil attribute plot
# @param df Data frame containing data for a field
# @param attribute Soil attribute
# @return Custom soil attribute plot
create_plot <- function(df, attribute) {
  map <- get_map(location = c(median(df$longitude), 
                              median(df$latitude)),
                 maptype = "satellite",
                 source = "google",
                 crop = FALSE,
                 zoom = 15)
  plot <- ggmap(map) + geom_point(aes_string(x = "longitude", 
                                             y = "latitude",
                                             color = attribute),
                                  size = 5,
                                  data = df)
  plot <- plot + ggtitle(paste("Farmer", df$farmer, "/ Field", df$field))
  plot <- plot + scale_color_gradient(low = "darkorange", high = "darkorchid4")
  return(plot)
}
ph_plot <- dlply(rice, "field", create_plot, attribute = "ph")
ph_plots <- do.call(arrangeGrob, ph_plot)

6a010534b1db25970b01bb07d2c8c1970d-800wi