Using R to find the “remotest” place in the world

For each of us, there exists a place on Earth which is the “remotest” place for us. That is, given all the places we have visited, which point on the globe is farthest away from everywhere we have ever been.

How do we find that place? We can use R’s geospatial libraries for that.

Here’s the core idea:

Step 1: We list all the places we’ve visited  (Lat, Long of those points.) Let’s call that set V.
Step 2: For any random point P on the globe, we can calculate its “remoteness” to us. That is the distance from P to the closest point in our visited set. For each v in V, we calculate its great circle distance to P and take the minimum.
Step 3: We do this for a large set of P that spans in the globe. (In say 1 degree intervals.)
Step 4: The point which has the greatest “remoteness” (max) is the point we are looking for.
This is sometimes called a “maximin” problem, the maximum of all minimum distances.

The idea is illustrated in this picture:

Screen Shot 2013-04-15 at 9.19.32 PM

Step 1: Write a few utility functions

Calculating Great Circle Miles between two points on Earth. For this I relied on Mario Pineda-Krch’s excellent post in r-bloggers. I just took the 3 functions I needed for this exercise.

# Code from
# Convert degrees to radians
deg2rad <- function(deg) return(deg*pi/180)

# Calculates the geodesic distance between two points specified by radian latitude/longitude using the
# Spherical Law of Cosines (slc)
gcd.slc <- function(long1, lat1, long2, lat2) {
  #  R <- 6371 # Earth mean radius [km]
  R <- 3961 # Earth mean radius [miles]  
  d <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R
  return(d) # Distance in km

gcd <- function(long1, lat1, long2, lat2) {
  # Convert degrees to radians
  long1 <- deg2rad(long1)
  lat1 <- deg2rad(lat1)
  long2 <- deg2rad(long2)
  lat2 <- deg2rad(lat2)     
  return(gcd.slc(long1, lat1, long2, lat2)  )  

Step 2: Create a vector of places visited

I found a list of 250 world cities, with their lat/long values, and then created a copy of that with a 0/1 marking the places I had visited. This gist has a working copy if you want to try it for the cities you have visited.

Step 3: Generate a list of points on the globe

Latitude: From 90N to 90S
Longitude: From 0 to 180E to 180E (via 0)

#Create a datafram of Long Lats that spans the Globe
getGlobespan <- function() {
  latsteps = seq(-90,90, by=10) #in deg
  lonsteps = seq(-180,180, by=10)
  df <- NULL
  for (lat in latsteps) {
    for (lon in lonsteps) {  
      #print(paste(lon,"deg", lat))
      df <- rbind(df, c(lon, lat))
  df <-
  names(df) <- c("Long", "Lat")

Step 4: Find the remotest point

Now, it is just a matter of calculating the “remoteness” of each point on the globe, and then taking its maximum.

unif.pts <- getGlobespan()
nearlist <- NULL
pt.closeness <- NULL
remote.dist <- 0
for (p in 1:nrow(unif.pts) ) {
  ptlong <- unif.pts[p,1]
  ptlat <- unif.pts[p,2]
  nearlist <- getMinDistanceFromOnePointToAllVisited(ptlong, ptlat, vdf)
  pt.closeness[p] <- nearlist[1] <- nearlist[1]
  if(remote.dist < { 
    remote.dist <-
    remote.long <- ptlong  <- ptlat
    closest.vlong  <- nearlist[2]
    closest.vlat   <- nearlist[3]
    print(c("remotest", ptlong,ptlat, round(remote.dist,2), closest.vlong, closest.vlat ))

Step 5: Plotting

Plotting these steps actually makes it very clear to see what’s going on.

First, on the world map, let’s plot all points – at 10 degree separations: This will be our proxy for “all” points on the globe.

mp <- NULL	 	
mapWorld <- borders("world", colour="gray70", fill="gray70")
title <- ggtitle("All points in the world, at 10 degree separation")	
uniform.pts <- geom_point(aes(x=gx, y=gy) ,color="black")
mp <- ggplot() +   mapWorld + title + uniform.pts	 


Similarly, we can plot all the visited points (in blue)

mp <- NULL	 
mapWorld <- borders("world", colour="gray70", fill="gray70")
title <- ggtitle("Plotting Places Visited")
visited.pts <- geom_point(aes(x=visitx, y=visity) ,color="blue")	 
mp <- ggplot() +   mapWorld + title + visited.pts


By adding a couple of lines, we can highlight the remotest point.

mp <- mp+ geom_point(aes(x=remote.long, , color="red", size=3 )
mp <- mp + geom_segment(aes(x=remote.long,, xend=closest.vlong, yend=closest.vlat))


Color-coding every point by its “remoteness”

One last thing: We can actually take *every* point on the globe and color-code by its remoteness. This serves to highlight the regions that haven’t been visited at all. The key line is the one that has “color = pt.closeness”

mp <- NULL
mapWorld <- borders("world", colour="gray70", fill="gray70")
title <- ggtitle("Every point on Earth color coded by its 'remoteness'")
mp <- ggplot() +   mapWorld + title

#Layer other points on top
gx <- unif.pts[,1]
gy <- unif.pts[,2]

#user closenss to color and determine square size. 15=square
mp <- mp+ geom_point(aes(x=gx, y=gy, color=pt.closeness, size=pt.closeness), shape=15) 
mp <- mp + scale_size(c(1:8))
mp <- mp + scale_color_gradient(high="red", low="green")