2015-12-03 66 views
2

我想知道是否有人用ggmap繪製經緯度點附近的圓半徑?例如,我想繪製一個給定的點,然後在這個點的2,500英尺的半徑內繪製陰影。我的腦海裏有一個想法,用大圓圈的公式來做這件事,但我想先看看這裏。ggmap創建的圖上緯度/長點附近的圓周半徑

回答

1

我一直在構建一個簡單的黑客,並不是我想象中的那個圈子,但它現在可以工作。

library(ggmap) 
##__________________________________________________________________ 
### earth.dist I found on r-bloggers. I believe it now belongs to the fossil package. 
earth.dist <- function (lat1,long1,lat2, long2) 
{ 
    rad <- pi/180 
    a1 <- lat1 * rad 
    a2 <- long1 * rad 
    b1 <- lat2 * rad 
    b2 <- long2 * rad 
    dlon <- b2 - a2 
    dlat <- b1 - a1 
    a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2 
    c <- 2 * atan2(sqrt(a), sqrt(1 - a)) 
    R <- 6378.145 
    d <- R * c 
    # this I changed to return feet instead of KM 
    return(d* 3280.8) 
} 

##__________________________________________________________________ 
## Function to output polygon to map 
BoxGon <- function(Lat,Long,feet){ 

    for(i in 1:1000){ 
    point = Long - i/1000 
    Dist <- earth.dist(Lat,Long, Lat,point) 
    if(Dist > feet){ 
        West <- cbind(Lat,point) 
        break} 
    } 

    for(i in 1:1000){ 
    point = Long + i/1000 
    Dist <- earth.dist(Lat,Long, Lat,point) 
    if(Dist > feet){ 
        East <- cbind(Lat,point) 
        break} 
    } 

    for(i in 1:1000){ 
    point = Lat + i/1000 
    Dist <- earth.dist(Lat,Long, point,Long) 
    if(Dist > feet){ 
        North <- cbind(point,Long) 
        break} 
    } 

    for(i in 1:1000){ 
    point = Lat - i/1000 
    Dist <- earth.dist(Lat,Long, point,Long) 
    if(Dist > feet){ 
        South <- cbind(point,Long) 
        break} 
    } 

    return(rbind(West,North,East,South,West)) 
} 

##__________________________________________________________________ 
df = BoxGon(37.295844, -121.898057,5000) 
df = as.data.frame(df) 
colnames(df) <- c('Latitude','Longitude') 

map <- get_map(location = 'san,jose', zoom = 12) 
map <- ggmap(map) 

# we select - 1 because once we map in pairs. IE once we have the last record there is nothing for that record to map to 
for(i in 1:nrow(df)-1){ 
latlon <- head(df,2) 
map <- map + geom_polygon(data=latlon,aes(x=Longitude,y=Latitude),alpha=0.1,size = 1,colour="green",fill="green") 
df = df[-1,] 
print(i) 
} 
map <- map + geom_polygon(aes(x=-121.898057,y=37.295844),alpha=0.1,size = 6,colour='Purple',fill='Purple') 
相關問題