Drawing a Circle with a Radius of a Defined Distance in a Map

nathanbeagle picture nathanbeagle · Apr 14, 2014 · Viewed 8.3k times · Source

I am able to plot a map and caption a specific point:

library(maps)
map("state")
text(-80.83,35.19,"Charlotte",cex=.6)

I can also plot a circle centered around that point:

symbols(-80.83,35.19,circles=2, add=TRUE)

However, I would like to control the size of the circle. In particular, I want to draw a circle with a radius of 100 mile around multiple locations contained in a data.frame, matrix, or list.

Answer

Gregoire Vincke picture Gregoire Vincke · Mar 18, 2015

The proposition of Gary is well adapted for plane maps, but can not be applied to maps generated by the "maps" package because it does not take care of the projection used to draw the map. Applied strictly as this, it results in drawing of an elipse (see below), because the unit of the circle radius is degrees, but not kilometers or miles. But degrees in latitude and longitude do not correspond to the same physical distance. To draw arround a point a circle, or something that is close to a circle, whose radius is a contant distance in miles or kilometers, you need to compute the corrected coordinates regarding to the projection. Taking your function and adapting it regarding to Chris Veness explanations on http://www.movable-type.co.uk, your function became :

library(maps)
library(mapdata)#For the worldHires database
library(mapproj)#For the mapproject function
plotElipse <- function(x, y, r) {#Gary's function ;-)
   angles <- seq(0,2*pi,length.out=360)
   lines(r*cos(angles)+x,r*sin(angles)+y)
}
plotCircle <- function(LonDec, LatDec, Km) {#Corrected function
    #LatDec = latitude in decimal degrees of the center of the circle
    #LonDec = longitude in decimal degrees
    #Km = radius of the circle in kilometers
    ER <- 6371 #Mean Earth radius in kilometers. Change this to 3959 and you will have your function working in miles.
    AngDeg <- seq(1:360) #angles in degrees 
    Lat1Rad <- LatDec*(pi/180)#Latitude of the center of the circle in radians
    Lon1Rad <- LonDec*(pi/180)#Longitude of the center of the circle in radians
    AngRad <- AngDeg*(pi/180)#angles in radians
    Lat2Rad <-asin(sin(Lat1Rad)*cos(Km/ER)+cos(Lat1Rad)*sin(Km/ER)*cos(AngRad)) #Latitude of each point of the circle rearding to angle in radians
    Lon2Rad <- Lon1Rad+atan2(sin(AngRad)*sin(Km/ER)*cos(Lat1Rad),cos(Km/ER)-sin(Lat1Rad)*sin(Lat2Rad))#Longitude of each point of the circle rearding to angle in radians
    Lat2Deg <- Lat2Rad*(180/pi)#Latitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
    Lon2Deg <- Lon2Rad*(180/pi)#Longitude of each point of the circle rearding to angle in degrees (conversion of radians to degrees deg = rad*(180/pi) )
    polygon(Lon2Deg,Lat2Deg,lty=2)
}
map("worldHires", region="belgium")#draw a map of Belgium (yes i am Belgian ;-)
bruxelles <- mapproject(4.330,50.830)#coordinates of Bruxelles
points(bruxelles,pch=20,col='blue',cex=2)#draw a blue dot for Bruxelles
plotCircle(4.330,50.830,50)#Plot a dashed circle of 50 km arround Bruxelles 
plotElipse(4.330,50.830,0.5)#Tries to plot a plain circle of 50 km arround Bruxelles, but drawn an ellipse

Result in image

(sorry my "reputation" do not allow me to post images ;-). Edit: Added your image.

I hope this helps. Grégoire