Equidistant points on a map

October 17, 2013
By

[This article was first published on Freakonometrics » R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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))
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))
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))
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))
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.

var vglnk = { key: '949efb41171ac6ec1bf7f206d57e90b8' };

(function(d, t) {
var s = d.createElement(t); s.type = 'text/javascript'; s.async = true;
s.src = '//cdn.viglink.com/api/vglnk.js';
var r = d.getElementsByTagName(t)[0]; r.parentNode.insertBefore(s, r);
}(document, 'script'));

Related
ShareTweet

To leave a comment for the author, please follow the link and comment on their blog:  Freakonometrics » R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or  here if you don't.

If you got this far, why not subscribe for updates from the site?  Choose your flavor: e-mail, twitter, RSS, or facebook...
```

Comments are closed.

Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)