I am trying to overlay two shapefiles on a raster grid. The raster is a digital
elevation model. The first shapefile is a polygon, which shows the outline of a
country, and the second shapefile is a polyline, which represents rivers. I am
currently using levelplot to plot the raster and spplot to overlay the polyline
and polygon. The polyline is coloured according to the attribute table using
the STATUS column - this plots "perennial" status as blue and "ephemeral"
status as orange. My code is shown below. This produces a plot with a colorkey
for the raster, but I can't work out how to add a key for the polyline.
Apologies if this is very obvious but I can't seem to find an answer on any
existing threads.
Any help would be much appreciated.
The code is as follows:
# Read rasterdata_location_DEM1 <-
"E:/AGA_ProcessDatasets/DEM_Data/gt30e020n40.tif"iDEM <-
raster(data_location_DEM1)
# Read shapefilesioutlineCountry <- readOGR(dsn = "E:/AGA_ProcessDatasets",
layer = "UNAfricanCountries_Final_Buffer100m")proj4string(ioutlineCountry) <-
CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
irivers <- readOGR(dsn = "E:/AGA_ProcessDatasets/Hydrology_Data", layer =
"Rivers_Africa_37333_ForPlotting_Split")proj4string(irivers) <-
CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
# define colours for riversrivercolours <- c("orange","blue")
# define colours for DEMcols <- gray.colors(128,start=1,end=0)
# define extents for plottingplotmin <- cellStats(iDEM,stat='min')plotmax <-
cellStats(iDEM,stat='max')breaks <- (plotmax - plotmin)/128legendbreaks <-
round((plotmax - plotmin)/5,0)scale <- seq(plotmin, plotmax, by=breaks)xmin <-
iDEM@extent@xminxmax <- iDEM@extent@xmaxymin <- iDEM@extent@yminymax <-
iDEM@extent@ymax
# define output png filefilenamePre =
paste(iCountryName,"_Hydrology2.png",sep="")png(filename="Hydrology.png",width=1000,height=486,units="px")
# plot
dataa<-levelplot(iDEM,margin=FALSE,xlab="",ylab="",at=scale,col.regions=cols,
colorkey=list(space="right",width=0.75,height=0.75)) +
spplot(irivers,zcol="STATUS",col.regions=rivercolours) +
spplot(ioutlineCountry,zcol="NAME_0",fill="transparent",col="black") # add
title to legend a$legend$right <- list(fun =
mergedTrellisLegendGrob(a$legend$right,list(fun = textGrob,args =
list("maSL",rot = 0)),vertical = FALSE))
# printprint(a,split=c(1,1,2,1),newpage=FALSE)
dev.off()
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
