Overlaying shapefile over multiple raster plots in R

Question

I have multiple NetCDF data and I need to plot only one component of each data. I have achieved this in plotting all the data. But now I am unable to overlay my shapefile over the plots created. Is there any way to overlay the shapefile over the maps I have plotted using levelplot? I tried it with spplot but still, I am not getting the shapefile over my maps. I also tried using sp.polygons but can still not make it. Here is a snippet of my code and the generated plot.

library(ncdf4)
library(raster)
library(rgdal)
library(rasterVis)
library(RColorBrewer)
library(colorRamps)
path<- setwd('D:/DATA/2015/MIXING_RATIO/')

# Set the path to the shapefile
shapefile_path <- "C:/Users/IITM/Downloads/India_State/India_State/India_State.shp"

# Load the shapefile
shapefile <- readOGR(shapefile_path)
# Get a list of all netCDF files in the directory
nc_files <- list.files(path, pattern = ".nc", full.names = TRUE)
raster_list <- list()
# Loop through each netCDF file
for (nc_file in nc_files) {
  # Read the netCDF data
  data <- nc_open(nc_file)
  
  # Extract the necessary variables
  lon <- ncvar_get(data, "lon")
  lat <- ncvar_get(data, "lat")
  dust_MR <- ncvar_get(data, "DU001")
  
  # Plot the first level slice of the data
  dust1.slice1 <- dust_MR[, , 60]
  
  # Create a raster object
  r <- raster(t(dust1.slice1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat), crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
  
  # Flip the raster vertically
  r <- flip(r, direction = 'y')
  
  # Crop the raster based on the shapefile extent
  r_crop <- crop(r, extent(shapefile))
  
  # Mask the raster based on the shapefile geometry
  r_mask <- mask(x = r_crop, mask = shapefile)
  
  # Add the raster to the list
  raster_list[[nc_file]] <- r_mask
  
  # Close the netCDF file
  nc_close(data)
}


# Stack the rasters
raster_stack <- stack(raster_list)
crs(raster_stack)
DUST<- stack(mget(rep("raster_stack")))
names(DUST)
raster_names<- gsub('D..DATA.2015.MIXING_RATIO.MERRA2_400.inst3_3d_aer_Nv.|.SUB.nc'," ",names(DUST))
raster_names
raster_dates <- as.Date(trimws(raster_names), format = "%Y%m%d")
raster_dates
txtdates <- format(raster_dates, "%d%b%Y")

plot(raster_stack, nc = 4)
# Create a single plot using levelplot
p<-levelplot(raster_stack,names.attr=txtdates,
          margin = FALSE,colorkey = list(space = "right"), main = "Dust Mixing Ratio (bin 1)(Kg/Kg)")
p + layer(sp.polygons(polygon, lwd=0.8, col='darkgray'))

Thank you for any suggestions.

score 1 · Accepted Answer · answered May 23 '23 at 18:53

1

This is another version of the answer by Robert Hijmans using rasterVis. The code includes a function, llines.SpatVector, that will be included in rasterVis after next update of the lattice package (more information here).

library(terra)
library(latticeExtra)
library(rasterVis)

llines.SpatVector <- function(x, ...) {
  xy <- crds(x, list=TRUE)
  names(xy) <- c("x", "y")
  lattice::llines(xy, ...)
}

r <- rast(system.file("ex/elev.tif", package="terra")) * 1:2
v <- vect(system.file("ex/lux.shp", package="terra"))
lns <- as.lines(v)

levelplot(r) + layer(llines(lns))

answered May 23 '23 at 18:53

Oscar Perpiñán

4,491
17
28

Thank you @Oscar it worked. But a problem is my plots are getting titled as you can see in the attached image file. How to tackle that? – shravani banerjee May 25 '23 at 05:30
@shravanibanerjee titled or tilted? Are you talking about [this question](https://stackoverflow.com/questions/76329168/lattice-subplots-appear-tilted)? – Oscar Perpiñán May 25 '23 at 09:18
sorry for the spell error. its tilted. – shravani banerjee May 25 '23 at 10:10

Robert Hijmans · Answer 2 · 2023-05-23T14:43:30.663

0

Here is how you can do that.

Example data

library(terra)
r <- rast(system.file("ex/elev.tif", package="terra")) * 1:2
v <- vect(system.file("ex/lux.shp", package="terra"))

Solution

plot(r, fun=\() lines(v))

Or if you want a single shared legend:

panel(r, fun=\() lines(v))

You can also use ggplot2 via the "tidyterra" package.

I should also note that the way you open the netCDF files is cumbersome and that the extent you set is not correct.

You should be able to simply do

nc <- terra::rast(nc_file)

(or the equivalent with the old "raster" package if you are still using that). And your entire script can probably be reduced to

library(terra)
v <- vect(shapefile_path)
nc_files <- list.files(path, pattern = ".nc", full.names = TRUE)
r <- rast(nc_files)
r <- r[[60 * (1:length(nc_files))]]
rcm <- crop(r, v, mask=TRUE)
panel(rcm, fun=\() lines(v))

edited May 23 '23 at 14:43

answered May 23 '23 at 14:32

Robert Hijmans

40,301
4
55
63

I need to select one variable and one dimension from my NC file and plot it. So using this nc <- terra::rast(nc_file), is not possible. – shravani banerjee May 24 '23 at 06:52
I think it is possible, as you can specify the variable of interest with `rast`. – Robert Hijmans May 24 '23 at 16:25
And then you would get the correct extent as well. – Robert Hijmans May 25 '23 at 09:49
Can you tell me how to select the variable of interest as well as the pressure dimension using rast? – shravani banerjee May 25 '23 at 12:10
I do not have your file so I do not know that. You can ask a new question with an example file. – Robert Hijmans May 26 '23 at 07:27

Overlaying shapefile over multiple raster plots in R

2 Answers2