3

I am trying overlay the shapefile of south asia on top of multiple raster plots using the code as below: 'a' is a multilayered raster file. Here is the link to the data (917KB size) Test_Data

ras <- list.files("/filepath/", pattern = "\\.tif$", full=TRUE)
s <- stack(ras)
south_asia  <- readOGR('/filepath/south_asia.shp')  #to import shapefile
cropped <- crop(x = s, y = extent(south_asia))          #crop raster
plot(cropped)
plot(south_asia, add=TRUE)

This code gives me one shapefile and multiple raster plots. How do i overlay the shapefile on top of the rasters? enter image description here Any help will be appreciated. P.S: They are in same CRS Thanks

Dipu
  • 123
  • 2
  • 14

2 Answers2

3

First, please try to give a reproducible example instead of a link to download files externally!

If you manually construct the plot (e.g. using par with base plotting) you can get your desired behaviour:

library(raster)

## testdata

# shapefile
shp <- getData(country='IND', level=1)

# raster 
r <- getData('alt', country='IND', mask=TRUE)

# create 4 layer rasterstack
rs <- stack(r,r,r,r)

## finally plot

# 2 rows, 2 cols
par(mfrow=c(2,2))


# loop layers
for (ii in 1:nlayers(rs)){

  plot(subset(rs,ii), main=names(rs)[ii])
  plot(shp, add=T)


}

Edit:

Use plot(subset(rs,ii), main=names(rs)[ii]) in the loop to plot the respective layer.

The result:

enter image description here

Val
  • 6,585
  • 5
  • 22
  • 52
  • Though your results is what i wanted. Still i am getting the same figure as above (mine). Please try the reproducible example as i edited the codes. – Dipu Jan 10 '20 at 08:50
  • The key in my approach is to use `par(mfrow=c(2,2))` to set up the plot and the visualize each layer independently instead of plotting the entire rasterstack at once. Now I see that I've included an error, as I plot the same raster over and over within the loop - instead it should be the layer `ii` in each iteration. I'll update – Val Jan 10 '20 at 08:52
  • ```plot(subset(rs,ii), main=names(rs)[ii])``` is doing the job but i get one plot on the screen. The last plot is the last layer on the screen. I want all the plots in one figure and not individual figures. Kindly update – Dipu Jan 10 '20 at 12:19
  • @Dipu did you run `par(mfrow=c(2,2))` prior? – Val Jan 10 '20 at 12:25
  • Thank you very much..i forgot to run that code Yes...it worked.. i was wondering if you can also suggest the code for one color bar for multiple plots and the color bar to be positioned below the figure with legend names with unites like Temperature (^oC) – Dipu Jan 10 '20 at 12:48
  • 1
    Things start to get a bit cumbersome using base plotting. You can do it, but maybe you want to check out [rasterVis](https://oscarperpinan.github.io/rastervis/) – Val Jan 10 '20 at 13:15
1

I would suggest using the argument addfun in the plot option for raster stack/brick For example:

# Function to add shapefile "shp_file" on each raster plot
add_shp=function(){plot(shp_file, bg="transparent", add=TRUE)}

#Plot selected raster of a raster stack, rs
plot(rs[[c(1:5)]],col=tim.colors(5),addfun=add_shp)

This will add shapefile shp_file to each of the 5 plots from the rasterstack rs.

Vinit Sehgal
  • 11
  • 2