title: How to efficiently calculate distance between GPS points in one dataset and GPS points in another data set using data.table
body: <p>I am facing a coding (optimization) problem in R. I have a long data set with GPS coordinates (lon, lat, timestamp) and for every row I need to check whether the location is near a bus stop. I have a .csv file with all the bus stops (in the Netherlands). The GPS coordinates file is millions of entries long, but could be split if necessary. The bus stop dataset is around 5500 entries long.<br>
  Using the code and tips given on, inter alia, these pages:</p>
  
  <p>1) <a href="https://stackoverflow.com/questions/36817423/how-to-efficiently-calculate-distance-between-pair-of-coordinates-using-data-tab">How to efficiently calculate distance between pair of coordinates using data.table :=</a></p>
  
  <p>2) <a href="https://stackoverflow.com/questions/29430303/using-a-simple-for-loop-on-spatial-data">Using a simple for loop on spatial data</a></p>
  
  <p>3) <a href="https://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula">Calculate distance between two latitude-longitude points? (Haversine formula)</a></p>
  
  <p>4) <a href="https://stackoverflow.com/questions/52910965/fastest-way-to-determine-country-from-millions-of-gps-coordinates-r">Fastest way to determine COUNTRY from millions of GPS coordinates [R]</a></p>
  
  <p>I was able to construct a code that works, but is (too) slow. I was wondering if someone can help me with a faster data.table() implementation or can point out where the bottle neck in my code is? Is it the spDistsN1() function,  or maybe the apply and melt() functions combination? I am most comfortable in R, but open to other software (as long as it is open source).</p>
  
  <p>Due to privacy concerns I cannot upload the full dataset, but this is a (small) reproducible example that is not too different from how the real data looks.</p>
  
  <pre><code># packages:
  library(data.table)
  library(tidyverse)
  library(sp)
  
  
  # create GPS data
  number_of_GPS_coordinates &lt;- 20000
  set.seed(1)
  gpsdata&lt;-as.data.frame(cbind(id=1:number_of_GPS_coordinates, 
                               lat=runif(number_of_GPS_coordinates,50.5,53.5), 
                               lon=runif(number_of_GPS_coordinates,4,7)))
  
  # create some busstop data. In this case only 2000 bus stops
  set.seed(1)
  number_of_bus_stops &lt;- 2000
  stop&lt;-as.data.frame(gpsdata[sample(nrow(gpsdata), number_of_bus_stops), -1]) # of course do not keep id variable
  stop$lat&lt;-stop$lat+rnorm(number_of_bus_stops,0,.0005)
  stop$lon&lt;-stop$lon+rnorm(number_of_bus_stops,0,.0005)
  busdata.data&lt;-cbind(stop, name=replicate(number_of_bus_stops, paste(sample(LETTERS, 15, replace=TRUE), collapse="")))
  
  names(busdata.data) &lt;- c("latitude_bustops",  "longitude_bustops", "name")
  </code></pre>
  
  <p>Download the real bus stop data if you want, kind of hard to reproduce a random sample of this.</p>
  
  <pre><code>#temp &lt;- tempfile()
  #download.file("http://data.openov.nl/haltes/stops.csv.gz", temp) #1.7MB
  #gzfile(temp, 'rt')
  #busstopdata &lt;- read.csv(temp, stringsAsFactors = FALSE)
  #unlink(temp)
  #bus_stops &lt;- fread("bus_stops.csv")
  #busdata.data &lt;- busstopdata %&gt;%
  #  mutate(latitude_bustops = latitude)%&gt;%
  #  mutate(longitude_bustops = longitude)%&gt;%
  #  dplyr::select(name, latitude_bustops,  longitude_bustops)
  </code></pre>
  
  <p>Code I use now to calculate distances. It works but it is pretty slow</p>
  
  <pre><code>countDataPoints3 &lt;- function(p) {
    distances &lt;- spDistsN1(data.matrix(gpsdata[,c("lon","lat")]), 
                           p,
                           longlat=TRUE) # in km
    return(which(distances &lt;= .2)) # distance is now set to 200 meters
  }
  
  
  # code to check per data point if a bus stop is near and save this per bus stop in a list entry
  datapoints.by.bustation       &lt;- apply(data.matrix(busdata.data[,c("longitude_bustops","latitude_bustops")]), 1, countDataPoints3)
  
  
  # rename list entries
  names(datapoints.by.bustation) &lt;- busdata.data$name
  
  # melt list into one big data.frame
  long.data.frame.busstops       &lt;- melt(datapoints.by.bustation)
  
  # now switch to data.table grammar to speed up process
  # set data.table
  setDT(gpsdata)
  gpsdata[, rowID := 1:nrow(gpsdata)]
  setkey(gpsdata, key = "rowID")
  setDT(long.data.frame.busstops)
  
  # merge the data, and filter non-unique entries 
  setkey(long.data.frame.busstops, key = "value")
  GPS.joined        &lt;- merge(x = gpsdata, y = long.data.frame.busstops, by.x= "rowID", by.y= "value", all.x=TRUE)
  GPS.joined.unique &lt;- unique(GPS.joined, by="id") # mak
  
  # this last part of the code is needed to make sure that if there are more than 1 bus stop nearby it puts these bus stop in a list
  # instead of adding row and making the final data.frame longer than the original one
  GPS.joined.unique2 &lt;- setDT(GPS.joined.unique)[order(id, L1), list(L1=list(L1)), by=id]
  GPS.joined.unique2[, nearby := TRUE][is.na(L1), nearby := FALSE] # add a dummy to check if any bus stop is nearby.
  
  # makes sense:
  as.tibble(GPS.joined.unique2) %&gt;%
    summarize(sum = sum(nearby)) 
  </code></pre>