Replicating QGIS color ramp style in R?

Question

I want to use rgdal to plot a shapefile graduating the colour of features according to the values of one field (density). I couldn't find online detailed indications on how to do it and the package documentation doesn't explain the options for plot().

This is what I composed so far

library(rgdal)  

shapefile = "shapefile.shp"
map <- readOGR(dsn = shapefile, layer = "shapefile", verbose=FALSE)

plot(map, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), border=NA, add=TRUE)

This is the style in QGIS I am trying to replicate:

Answer 1

Have a look at spplot. Here is a short example using population density data from Haiti. I like to use Colorbrewer for creating color palettes suitable for maps, but there surely exist a lot of other possibilities to create sequential color scales in R.

# Required packages
library(rgdal)
library(RColorBrewer)

# Import shapefile
shp <- readOGR(dsn = "/path/to/data", layer = "Haiti_ADM3_stats")

# Color palette
my.colors <- brewer.pal(9, "Reds")

# Plot population density
spplot(shp, "POP_DENS", col.regions = my.colors, cuts = 8, 
       scales = list(draw = T))

Update:

Indeed, spplot takes some time to generate the graph. If you are looking for a more convenient way to plot your data as concerns speed, I would suggest to have a look at the ggplot package. Note that the following code is strongly influenced by this wordpress post on speeding up map plotting in R.

# Required package
library(ggplot2)

# Extract polygon corners and merge with shapefile data
shp@data$id <- rownames(shp@data)
shp.ff <- fortify(shp)
shp.df <- merge(shp@data, shp.ff, by = "id", all.y = T)

# Plot population density
ggplot() + 
  geom_polygon(data = shp.df, aes(x = long, y = lat, group = group, 
                                  fill = POP_DENS), color = "black") + 
  scale_fill_gradient(name = "Population density \n (inhabitants / km²)", 
                      low = my.colors[1], high = my.colors[9]) + 
  labs(x = "Lon", y = "Lat") + 
  theme_bw()

Answer 2

It seems you are working with the same data you mentioned on a previous question of yours.

This is what I got with R package ggplot2.

library(raster)
library(rgdal)
library(ggplot2)
library(maptools)
library(plyr)
library(gpclib)

#import shapefile
italy_map <- readOGR(dsn = "C:\\...\\ITA_adm", layer = "ITA_adm1")

#prepare shapefile for ploting in ggplot2, according to this page:  https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles
italy_map@data$id = rownames(italy_map@data)
gpclibPermit() #avoid "Error: isTRUE(gpclibPermitStatus()) is not TRUE" in the next command
italy_map.points = fortify(italy_map,region="id")
italy_map.df = join(italy_map.points, italy_map@data, by="id")

#import raster data for Italy population density in 2000
pop_density_file_2000 ="C:\\...\\ita_gpwv3_pdens_wrk_25\\itadens\\itads00ag\\w001001.adf"
dens_2000 <- raster(pop_density_file_2000)

#prepare raster data to plot with ggplot2
raster_points = rasterToPoints(dens_2000)
raster_df = data.frame(raster_points)

cuts = c(10,1000,5000,10000,15000) # set vector for legend breaks
legend_title="Hab/km²" #name legend title

ggplot() + coord_equal() + theme_bw() + 
       geom_tile(data=raster_df,aes(x=x,y=y,fill=w001001)) + 
       scale_fill_gradient(legend_title, low="white", high="red", breaks=cuts) + 
       geom_path(data=italy_map.df,aes(long,lat,group=group),colour="black") + 
       xlab("Longitude") + ylab("Latitude")