[R] To map population of European countries from Eurostat

[Arnaud Michel]

Hello
I would like to map the population of the European countries in 2011.
I am using the spatial shapefiles of Europe published by EUROSTAT.
I applyed the script below of Markus Kainubut but I had a problem with 
the map.

Any ideas ?
Thanks for your help

#########################################
# Importing the data
library(SmarterPoland)
df <- getEurostatRaw(kod = "tps00001")

# rename 1ere variable
names(df) <- c("xx", 2002:2013)

# new variable df$unit et df$geo.time
df$unit <- lapply(strsplit(as.character(df$xx), ","), "[", 1)
df$geo.time <- lapply(strsplit(as.character(df$xx), ","), "[", 2)

# file df.l
df.l <- melt(data = df, id.vars = "geo.time", measure.vars = c("2002", 
"2003",
"2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012", 
"2013"))

df.l$geo.time <- unlist(df.l$geo.time)  # unlist the geo.time variable

# Load the shapefile GISCO, subset it to NUTS2-level
# and merge the Eurostat
# attribute data with it into single SpatialPolygonDataFrame

Filetodownload <- 
"http://epp.eurostat.ec.europa.eu/cache/GISCO/geodatafiles/NUTS_2010_60M_SH.zip"
download.file(Filetodownload, destfile = "NUTS_2010_60M_SH.zip")
rm(Filetodownload)

# Unzip NUTS_2010_60M_SH.zip to SpatialPolygonsDataFrame
unzip("NUTS_2010_60M_SH.zip")

library(rgdal)
# Lecture du fichier
map <- readOGR(dsn = "./NUTS_2010_60M_SH/data", layer = "NUTS_RG_60M_2010")
# OGR data source with driver: ESRI Shapefile
# Source: "./NUTS_2010_60M_SH/data", layer: "NUTS_RG_60M_2010"
# with 1920 features and 4 fields
# Feature type: wkbPolygon with 2 dimensions

names(map)
#[1] "NUTS_ID"    "STAT_LEVL_" "SHAPE_Leng" "SHAPE_Area"

# as the data is at NUTS2-level, we subset the
# spatialpolygondataframe
map_nuts2 <- subset(map, STAT_LEVL_ <= 2)

# dim show how many regions are in the spatialpolygondataframe
dim(map_nuts2)
## [1] 467   4

# dim show how many regions are in the data.frame
dim(df)
## [1] 43  15

################################################################
# Spatial dataframe has 467 rows and attribute data 43.  We need to make
# attribute data to have similar number of rows
NUTS_ID <- as.character(map_nuts2$NUTS_ID)
VarX <- rep("empty", 467)
dat <- data.frame(NUTS_ID, VarX)

# then we shall merge this with Eurostat data.frame
dat2 <- merge(dat, df, by.x = "NUTS_ID", by.y = "geo.time", all.x = TRUE)

## merge this manipulated attribute data with the spatialpolygondataframe
## there are still duplicates in the data, remove them
dat2$dup <- duplicated(dat2$NUTS_ID)
dat3 <- subset(dat2, dup == FALSE)

## rownames
row.names(dat3) <- dat3$NUTS_ID
row.names(map_nuts2) <- as.character(map_nuts2$NUTS_ID)

## order data
dat3 <- dat3[order(row.names(dat3)), ]
map_nuts2 <- map_nuts2[order(row.names(map_nuts2)), ]

## join
library(maptools)
shape <- spCbind(map_nuts2, dat3)

##########################################################################
# Munging the shapefile into data.frame and ready for ggplot-plotting

## fortify spatialpolygondataframe into data.frame
library(ggplot2)
library(rgeos)
shape$id <- rownames(shape at data)
map.points <- fortify(shape, region = "id")
map.df <- merge(map.points, shape, by = "id")

###############################################################
# Only year 2011
# As we want to plot map faceted by years from 2003 to 2011 we have to
# melt it into long format
map.df <- map.df[, -c(15:23,25,26)]
library(reshape2)

map.df.l <- melt(data = map.df, id.vars = c("id", "long", "lat", 
"group"), measure.vars = c("X2011"))

# year variable (variable) is class string and type X20xx.  Lets remove
# the X and convert it to numerical

library(stringr)
map.df.l$variable <- str_replace_all(map.df.l$variable, "X", "")
map.df.l$variable <- factor(map.df.l$variable)
map.df.l$variable <- 
as.numeric(levels(map.df.l$variable))[map.df.l$variable]

# And finally the plot using ggplot2
#Map shows proportion of materially deprived households at the NUTS2 level.
#Grey color indicates missing data.


library(ggplot2)

ggplot(map.df.l, aes(long,lat,group=group)) +
geom_polygon(aes(fill = value)) +
geom_polygon(data = map.df.l, aes(long,lat),
fill=NA, color = "white",
size=0.1) + # white borders
coord_map(project="orthographic", xlim=c(-22,34), ylim=c(35,70)) + # proj
labs(title = "Cartographie") +
theme_minimal()

-- 
Michel ARNAUD
Chargé de mission auprès du DRH
DGDRD-Drh - TA 174/04
Av Agropolis 34398 Montpellier cedex 5
tel : 04.67.61.75.38
fax : 04.67.61.57.87
port: 06.47.43.55.31