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


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:

enter image description here


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

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



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

# 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") + 


  • This is probably a good solution. But I noticed the ssplot function is extremely slow in composing the plot with 8000 features... – Francesco Nov 23 '13 at 22:10
  • As concerns speed, you're definitely right! Maybe the above update helps ;-) – fdetsch Nov 24 '13 at 7:04

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.


#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")

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.