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

April 16, 2013 Leave a comment

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:

## 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

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

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

## 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