At first I want to show how you can make simple maps with R. Furthermore I want to demonstrate one important function from the package “gdalUtils”, which contains wrappers for the Geospatial Data Abstraction Library (GDAL) Utilities. I  will rasterize a vector with the function gdal_rasterize, because  it is quite fast for big vector data sets.

library(ggplot2)
library(maptools)
library(rgdal)
library(raster)
library(plyr)
library(gdalUtils)
library(RColorBrewer)
library(rasterVis)

Load world data (package maptools) and plot it to get an over view.

data(wrld_simpl)
plot(wrld_simpl)

Subset data. Select country of Austria (aut_poly) & Austria with neighbors in the west (aut_neigh).

aut<-wrld_simpl[wrld_simpl@data$NAME=="Austria",]
aut_neigh<-wrld_simpl[wrld_simpl@data$SUBREGION==155,]

Ggplot only supports data.frames, that is the reason why we transform the spatial polygon class to a data frame.

aut_gg<-fortify(aut)
aut_neigh_gg<-fortify(aut_neigh)
p1<-ggplot()+geom_map(data=aut_neigh_gg,map=aut_neigh_gg,aes(x=long,y=lat,map_id=id),fill="white",color="black", size=0.25)+geom_map(data=aut_gg,map=aut_gg,aes(map_id=id,fill=id))
p1

Now we change the theme of the map and add labels. Labels need a vector to define the
position. Therefore we calculate the mean X- and Y-coordinate and save both in the data.frame “labs.out_neigh”.
p2 shows us a new map with a title and labels for each country. If you want a better position for the lables you can expand the function geom_text with further arguments: Some  arguments are angle, vjust, size….

labs.aut_neigh<-ddply(aut_neigh_gg,.(id),summarise, meanx = mean(range(long)),meany=mean(range(lat)))
p2<-p1+ggtitle("Part of Europe")+geom_text(data=labs.aut_neigh,aes(meanx,meany,label=id))
p3<-p2+theme_bw()+theme(legend.position = "right")+
scale_fill_discrete(name="Country-Id",l=50)
p3

I am often faced with the task to transform vector data (shape files) to a raster file. There are few ways to get a solution. I have prepared two examples.

Solution one :

extent(aut_neigh) # Delivers us the extent of our area of interest.
raster_aut_neigh_1<-raster(extent(aut_neigh)) # Creates an raster based on the extent
raster_aut_neigh_1 # Shows use the number of rows and columns. Used in the following row
raster_aut_neigh_1[] = runif(10*10) # Fills the rasters with randomly selected values. 
plot(raster_aut_neigh_1) # Plot raster
plot(aut_neigh,add=TRUE) # Add polygon to the plot

res(raster_aut_neigh_1)=c(0.5,0.5) # We increase the spatial resolution
raster_aut_neigh_1[] = runif(28*44) 
plot(raster_aut_neigh_1)
plot(aut_neigh,add=TRUE)

raster_aut_neigh_final<-rasterize(aut_neigh,raster_aut_neigh_1,field="ISO3")
proj4string(raster_aut_neigh_final)=aut_neigh@proj4string
plot(raster_aut_neigh_final)

Solution two:

Let’s assume that we already have a shape file. So we use the gdalUtils package and transform the shape file to a raster file. At first we create a shape file of our area of interest.

diroutput<-"C:/Users/YourName/Documents/Foo/World_map" #define your path
# Writes a shape file into your directory
writeOGR(aut_neigh,dsn=diroutput,layer="aut_neigh",driver="ESRI Shapefile") 
# Load shape file in R
shape_file<-readOGR(dsn=diroutput,"aut_neigh")

Now we transform the shape file  (shape_file) in to a raster file.

#find EPSG code for the Raster
EPSG <- make_EPSG() # Functions to create a data.frame which contains all EPSG-Codes
aut_neigh@proj4string@projargs # obtains the projection
# Looking after the EPSG-Codes: you will get one row with EPSG: 4326
EPSG[grepl("+proj=longlat", EPSG$prj4, fixed=TRUE) & grepl("+datum=WGS84", EPSG$prj4, fixed=TRUE), ]
#> 249 4326 # WGS 84 +proj=longlat +datum=WGS84 +no_defs

Finally we define two paths for the last calculation.  Please adapte my code to your system properties.

inputdir<-"C:/Users/Foo/Documents/Ba/World_map"
outputdir<-"C:/Users/Foo/Documents/Ba/World_map/gdal_shp_2_r.tif"

#tr= target layer resolution; a_srs=EPSG-Code; l= shape file layer name; 
#a = shape file attribute field to copy
r_id<-gdal_rasterize(inputdir, outputdir,tr=c(0.5,0.5),a_nodata=-9999,a_srs="epsg:4326",l="aut_neigh",a="UN",output_Raster = FALSE)
#load rasterized polygon
r_id2<-raster(outputdir)
# define raster field attribute name
names(r_id2)<-"ID"

Now we plot our final results.

# We take a look to the relationship between country and the raster cell value "UN"
sort_aut<-aut_neigh@data[order(aut_neigh@data$UN),]

With following commands we reclassify the numeric raster call values.
Instead of the numeric cell values the legend shows us categories, which relate to the numeric attribute “UN”

r_test <- ratify(r_id2) 
rat <- levels(r_test)[[1]]
rat$country <- c(as.character(sort_aut[sort_aut$UN%in%rat[,1],"ISO3"]))
levels(r_test)<-rat
my.color<-brewer.pal(7,"Dark2") # define your own color palette
levelplot(r_test,col.regions=my.color)