Equidistant points on a map

This morning, I had a comment on a recent post, regarding a graph I did upload on the blog, which was extracted from a paper now online (see http://hal.archives-ouvertes.fr/hal-00871883). Jo (from KUL, I guess I can share that piece of information) asked me

I was wondering whether you would want to share the R code for plotting figures 1 and 14? W.r.t. the former, the figure-in-figure is a nice touch; as to the latter, I am curious to know how you translated distance in km to the size parameters of the graph (par(“usr”)) for plotting the corresponding concentric circles (and the arrow indicating the radius) on top of your map.

At first, I thought I made a mistake in my plot. I mean, each time I have a question, I start to be suspicious, and I start to wonder if what I did was valid, or not. Here was the graph

Let’s make it clear: I do not draw circles here. So yes, I believe that what I did is valid. What I did is simple. First, I get the background map,

library(maps)
map("world",xlim=c(130,150),ylim=c(25,45),fill=TRUE,col="light green")

Then, I need some function to compute distance from coordinates. The functions I use are

deg2rad = function(deg) return(deg*pi/180)
DISTANCEDEG = function(long1, lat1, long2, lat2) {
R=6371; d=acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R
return(d) 
}

The center here will be Tokyo (東京),

X=139+45/60
Y=35+40/60

The idea now is simple: I generate a grid (here 501×501)

VX=seq(X-10,X+10,length=501)
VY=seq(Y-10,Y+10,length=501)
VtX=rep(VX,each=length(VY))
VtY=rep(VY,length(VX))
ZDeg=deg2rad(cbind(VtX,VtY))

I compute the distance from all those points to Tokyo, and check is the distance is larger or smaller than a given value,

L=500
D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L

If the distance is smaller than 500km, then I put a blue dot on the graph,

points(VtX[D1],VtY[D1],pch=19,cex=.2,col="light blue")

Then I use the same procedure for 250km (obviously, it is more convenient to start from larger and to go to smaller distances)

L=250
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtX[D],VtY[D],pch=19,cex=.2,col="light yellow")

Then, I did draw an arrow to ilustrate the largest distance

k=which.max(VtX[D1])
arrows((VtX[D1])[k],(VtY[D1])[k],X,Y,code=3,length=.1)
text(((VtX[D1])[k]+X)/2,Y+.35,"500 km")

And now, I have the graph.

Now, the point is that it should depend on the kind of projection we use, right? So here is a function that can be used for different kind of projections (some slight changes are necessary, since the map is now centered on some point, and we cannot use standard coordinates)

library(mapproj)
mapjapan = function(pr="conic",pm=45){
map("world","japan",fill=TRUE,col="light green",projection=pr, par=pm)
MP=mapproject(data.frame(x=X,y=Y),projection="")
Xp=MP$x
Yp=MP$y
VX=seq(X-10,X+10,length=501)
VY=seq(Y-10,Y+10,length=501)
VtX=rep(VX,each=length(VY))
VtY=rep(VY,length(VX))
MP=mapproject(data.frame(x=VtX,y=VtY),projection="")
VtXp=MP$x
VtYp=MP$y
ZDeg=deg2rad(cbind(VtX,VtY))
L=500
D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue")
L=250
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
L=100
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue")
L=50
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
points(Xp,Yp,pch=19,cex=.4,col="red")
map("world","japan",projection=pr, par=pm,add=TRUE)
}

The default function here produces a map based on a conic projection,

mapjapan()

But we can also use a Bonne projection (a pseudo-conic one, named after Rigobert Bonne)

mapjapan("bonne")

or Lagrange projection,

mappjapan("lagrange",NULL)

and (as a last one), Albers projections,

mapjapan("albers",c(30,40))

Of course, much more projections are possible !

We do not see much here, right ? So let us play with a larger country to visualize something. Like Canada. And the distance to, say, Winnipeg.

The first thing to do is to define the coordinates of Winnipeg,

X=-(97+08/60)
Y=(49+53/60)

Then, we slightly change our function

mapcanada = function(pr="conic",pm=45){
map("world","canada",fill=TRUE,col="light green",projection=pr, par=pm)
MP=mapproject(data.frame(x=X,y=Y),projection="")
Xp=MP$x
Yp=MP$y
VX=seq(X-30,X+30,length=501)
VY=seq(Y-30,Y+30,length=501)
VtX=rep(VX,each=length(VY))
VtY=rep(VY,length(VX))
MP=mapproject(data.frame(x=VtX,y=VtY),projection="")
VtXp=MP$x
VtYp=MP$y
ZDeg=deg2rad(cbind(VtX,VtY))
L=2000
D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue")
L=1000
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
L=500
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue")
L=200
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
points(Xp,Yp,pch=19,cex=.4,col="red")
map("world","canada",projection=pr, par=pm,add=TRUE)
}

Now, we can have some fun

mapcanada()

mapcanada("bonne",45)

mapcanada("albers",c(30,40))

mapcanada("lagrange",NULL)

Fun, isn’t it ? Changing the projection will change the shape of equidistant curves.