2017-12-11 65 views
5

bir ggmap üretilen rasterin üstüne bir geom_sf üretilen harita koymak için:Aşağıdaki kod çalıştı

library(ggplot2) 
library(ggmap) 
library(sf) 
nc <- st_read(system.file("shape/nc.shp", package = "sf")) 
str(nc) 

Classes ‘sf’ and 'data.frame': 100 obs. of 15 variables: 
$ AREA  : num 0.114 0.061 0.143 0.07 0.153 0.097 0.062 0.091 0.118 0.124 ... 
$ PERIMETER: num 1.44 1.23 1.63 2.97 2.21 ... 
$ CNTY_ : num 1825 1827 1828 1831 1832 ... 
$ CNTY_ID : num 1825 1827 1828 1831 1832 ... 
$ NAME  : Factor w/ 100 levels "Alamance","Alexander",..: 5 3 86 27 66 46 15 37 93 85 ... 
$ FIPS  : Factor w/ 100 levels "37001","37003",..: 5 3 86 27 66 46 15 37 93 85 ... 
$ FIPSNO : num 37009 37005 37171 37053 37131 ... 
$ CRESS_ID : int 5 3 86 27 66 46 15 37 93 85 ... 
$ BIR74 : num 1091 487 3188 508 1421 ... 
$ SID74 : num 1 0 5 1 9 7 0 0 4 1 ... 
$ NWBIR74 : num 10 10 208 123 1066 ... 
$ BIR79 : num 1364 542 3616 830 1606 ... 
$ SID79 : num 0 3 6 2 3 5 2 2 2 5 ... 
$ NWBIR79 : num 19 12 260 145 1197 ... 
$ geometry :sfc_MULTIPOLYGON of length 100; first list element: List of 1 
..$ :List of 1 
.. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ... 
    ..- attr(*, "class")= chr "XY" "MULTIPOLYGON" "sfg" 
    - attr(*, "sf_column")= chr "geometry" 
    - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ... 
    ..- attr(*, "names")= chr "AREA" "PERIMETER" "CNTY_" "CNTY_ID" ... 

map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google") 
a <- unlist(attr(map,"bb")[1, ]) 
bb <- st_bbox(nc) 
ggplot() + 
    annotation_raster(map, xmin = a[2], xmax = a[4], ymin = a[1], ymax = a[3]) + 
    xlim(c(bb[1], bb[3])) + ylim(c(bb[2], bb[4])) + 
    geom_sf(data = nc, aes(fill = AREA)) 

iki katman düzgün örtüşmeyen; coord_sf() ile projeksiyonu değiştirmeyi denedim, ancak başarılı olmadım.

herhangi bir öneri? teşekkürler

+0

sen str 'sonucunu gönderebilir miyim olduğu (nc) Sen Tahmin etmek gerekir koordinat başvuru sistemini belirtmek gerekir (ve doğru öylesine benziyor) (sorunuzu düzenleyin)? – Phil

cevap

2

Bunu yapmak için çizim yöntemini sf paketinden kullanabilirsiniz. 3857.

library(ggplot2) 
library(ggmap) 
library(sf) 

nc_shp <- st_read(system.file("shape/nc.shp", package = "sf")) 
nc_map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google") 

st_crs(nc_map) 
# Coordinate Reference System: NA 

# assume the coordinate refence system is 3857 
plot(st_transform(nc_shp, crs = 3857)[1], bgMap = nc_map) 

enter image description here

+0

çok teşekkür ederim. – user1453488

İlgili konular