Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Plotting a raster behind a shapefile

How can I plot a "raster" object behind a shapefile object? Both plot fine on their own but the points don't plot over the raster:

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)
like image 807
Benjamin Avatar asked Feb 14 '12 16:02

Benjamin


People also ask

How do you plot raster data?

Format Plot If we need to create multiple plots using the same color palette, we can create an R object ( myCol ) for the set of colors that we want to use. We can then quickly change the palette across all plots by simply modifying the myCol object. We can label the x- and y-axes of our plot too using xlab and ylab .

What does a raster plot show?

A raster plot is a simple method to visually examine the trial-by-trial variability of the respones. You can examine what features these responses have in common by averaging over all responses to create a peri-stimulus time histogram.


1 Answers

Overplotting raster plots with points, lines, and polygons should work just fine, as the following example shows.

My best guess would be that the Spatial* objects you are attempting to plot on top of the raster fall outside of the region being plotted. Have you checked that both the raster and Spatial* objects are in the same CRS, and (assuming they are) that the bounding boxes overlap? (i.e. try bbox(shp) and bbox(ras), and compare the results).

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

like image 132
Josh O'Brien Avatar answered Sep 24 '22 01:09

Josh O'Brien