2012-02-14 19 views
10

Jak wykreślić obiekt "rastrowy" za obiektem shapefile? Zarówno fabuła dobrze na własną rękę, ale punkty nie działki nad rastra:Wykreślanie rastra za shapefile

require(rgdal) 
require(maptools) 
require(raster) 

myproj = "+proj=utm +zone=12 +north +ellps=WGS84 +units=m" 
shp = readShapeSpatial(fn.shp, proj4string = CRS(myproj)) 
ras = raster(fn.tif) 

plot(ras) 
plot(shp, bg="transparent", add=TRUE) 

Odpowiedz

17

Overplotting działek rastrowe z punktów, linii i wielokątów powinien działać tak samo dobrze, jak w poniższym przykładzie.

Moim zdaniem domyślnie obiekty Spatial*, które próbujesz wydrukować na rastrze, wykraczają poza wykreślany obszar. Czy sprawdziłeś, czy obiekty raster i Spatial* są w tym samym CRS i (zakładając, że są), że pola obwiedniowe zachodzą na siebie? (tj. spróbuj bbox(shp) i bbox(ras) i porównaj wyniki).

library(rgdal) 
library(raster) 
# Create a raster 
ras <- raster(ncols=36, nrows=18) 
ras[] <- runif(ncell(ras)) 
# Create a SpatialPoints object 
shpPts <- spsample(Spatial(bbox=bbox(ras)), 20, type="random") 
# Create a SpatialPolygons object 
p1 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)) 
shpPolys <- SpatialPolygons(list(Polygons(list(Polygon(p1)), 1))) 

# Plot them, one layer after another 
plot(ras) 
plot(shpPts, pch=16, col="red", add=TRUE) 
plot(shpPolys, col="yellow", add=TRUE) 

enter image description here

+0

Dobra rozmowa, moja raster został nieprawidłowo położone. Dziękujemy za potwierdzenie, że działa to poprawnie. – Benjamin