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 http://www.r-bloggers.com/great-circle-distance-calculations-in-r/
# 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 <- as.data.frame(df)
  names(df) <- c("Long", "Lat")
  return(df)
}

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]
  closest.to.p <- nearlist[1]
  if(remote.dist < closest.to.p) { 
    remote.dist <- closest.to.p
    remote.long <- ptlong
    remote.lat  <- 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	 
mp

world

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
mp	

visited

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

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

remotest

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")
mp

color_coded_remoteness

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: