posted also on gis.stackexchange
I'm having some trouble plotting a RasterLayer and a SpatialPolygonsDataFrame in a level plot in R. Something is wrong with the projection, but I don't understand what.
Here you can find my reproducible example data and code:
library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(rasterVis)
library(maptools)
setwd("C:/...path_to/test2")
data<-read.csv("test.csv", header=TRUE)
#creating the raster from a data.frame and giving projection
raster<-rasterFromXYZ(data, crs="+proj=longlat")
raster_proj<-projectRaster(raster, crs="+proj=longlat +datum=WGS84
+ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
#loading shapes
nuts<-readShapePoly("NUTS_RG_WGS84.shp", proj4string=CRS("+proj=longlat"))
countries<-readShapePoly("countries.shp", proj4string=CRS("+proj=longlat"))
#subsetting the raster
raster_clip<- crop(raster_proj, countries, snap="near")
#plot
p.strip <- list(cex=1.5, lines=1, fontface='bold')
x.scale <- list(cex=1.5)
y.scale <- list(cex=1.5)
label <- list(labels=list(width=1, cex=1.5), height=0.95)
levelplot(raster_clip, par.settings = RdBuTheme, margin=FALSE,
at=seq(min(na.omit(values(raster_clip))), max(na.omit(values(raster_clip))), length.out=15),scales=list(x=x.scale, y=y.scale),
par.strip.text=p.strip, colorkey=label, xlab=list(label="Longitude", cex=2), ylab=list(label="Latitude", cex=2))+
layer(sp.polygons(nuts))+layer(sp.polygons(countries))
The error message is:
Error: Attempted to create layer with no geom.
But the raster and the polygons are in the same projection. I also tried to re-project the polygons using:
nuts <- spTransform(nuts, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"))
countries <- spTransform(countries, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"))
Any Idea?
If you don't want to or can't detach ggplot2 because you need it for something else, you can also call the layer function specifically from the latticeExtra package using:
levelplot(myRaster) +
latticeExtra::layer(sp.polygon(myPolygon))
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With