我有多个NetCDF数据，我只需要绘制每个数据的一个组件。我在绘制所有数据时做到了这一点。但现在我无法将我的shapefile覆盖在创建的图上。有没有办法将shapefile覆盖在我使用levelplot绘制的Map上？我尝试了它与spplot，但仍然，我没有得到shapefile超过我的Map。我也试过使用sp.多边形，但仍然不能使它。下面是我的代码片段和生成的图。

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'))

谢谢你的任何建议。