Gary 的命题非常适合平面地图,但不能应用于“maps”包生成的地图,因为它不考虑用于绘制地图的投影。严格如此应用,它会导致绘制椭圆(见下文),因为圆半径的单位是度,而不是公里或英里。但纬度和经度的度数并不对应于相同的物理距离。要围绕一个点绘制一个圆或接近圆的物体,其半径是以英里或公里为单位的恒定距离,您需要计算有关投影的校正坐标。采用您的功能并根据 Chris Veness 的解释进行调整http://www.movable-type.co.uk,你的函数变成了:
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
(抱歉,我的“声誉”不允许我发布图片;-)。编辑:添加了您的图像。
我希望这有帮助。
格雷瓜尔