Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

world map - map halves of countries to different colors using ggplot2

Tags:

plot

r

ggplot2

I'm looking for a ggplot2 solution for this question:

world map - map halves of countries to different colors

I reproduce the example from that question below, which is based on the question here (ggplot map with l).

library(rgdal)
library(ggplot2)
library(maptools)

# Data from http://thematicmapping.org/downloads/world_borders.php.
# Direct link: http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip
# Unpack and put the files in a dir 'data'

gpclibPermit()
world.map <- readOGR(dsn="data", layer="TM_WORLD_BORDERS_SIMPL-0.3")
world.ggmap <- fortify(world.map, region = "NAME")

n <- length(unique(world.ggmap$id))
df <- data.frame(id = unique(world.ggmap$id),
                 growth = 4*runif(n),
                 category = factor(sample(1:5, n, replace=T)))

## noise
df[c(sample(1:100,40)),c("growth", "category")] <- NA


ggplot(df, aes(map_id = id)) +
     geom_map(aes(fill = growth, color = category), map =world.ggmap) +
     expand_limits(x = world.ggmap$long, y = world.ggmap$lat) +
     scale_fill_gradient(low = "red", high = "blue", guide = "colorbar")
like image 913
Xu Wang Avatar asked Dec 21 '14 16:12

Xu Wang


2 Answers

You've got a couple options. It's pretty straightforward to get the polygons plotted, but you can't have two different fill scales. This solution uses point-shape annotations, but can be changed to scale the geom_point by color (or both color and shape). I think that's the best you're going to be able to do save overlaying two maps by hand in a separate program.

You'll also (probably) want to tweak the U.S. bounding box since the center is a bit off (a few of them are, actually, but that one is really obvious).

I also removed Antarctica. You can add that back in if you want, but it's wasted plot real estate IMO.

library(rgdal)
library(ggplot2)
library(maptools)
library(rgeos)
library(RColorBrewer)

world.map <- readOGR(dsn="/Users/bob/Desktop/TM_WORLD_BORDERS_SIMPL-0.3/", layer="TM_WORLD_BORDERS_SIMPL-0.3")

# Get centroids of countries
theCents <- coordinates(world.map)

# extract the polygons objects
pl <- slot(world.map, "polygons")

# Create square polygons that cover the east (left) half of each country's bbox
lpolys <- lapply(seq_along(pl), function(x) {
  lbox <- bbox(pl[[x]])
  lbox[1, 2] <- theCents[x, 1]
  Polygon(expand.grid(lbox[1,], lbox[2,])[c(1,3,4,2,1),])
})

# Slightly different data handling
wmRN <- row.names(world.map)

n <- nrow(world.map@data)
world.map@data[, c("growth", "category")] <- list(growth = 4*runif(n),
                                                  category = factor(sample(1:5, n, replace=TRUE)))

# Determine the intersection of each country with the respective "left polygon"
lPolys <- lapply(seq_along(lpolys), function(x) {
  curLPol <- SpatialPolygons(list(Polygons(lpolys[x], wmRN[x])),
                             proj4string=CRS(proj4string(world.map)))
  curPl <- SpatialPolygons(pl[x], proj4string=CRS(proj4string(world.map)))
  theInt <- gIntersection(curLPol, curPl, id = wmRN[x])
  theInt
})

# Create a SpatialPolygonDataFrame of the intersections
lSPDF <- SpatialPolygonsDataFrame(SpatialPolygons(unlist(lapply(lPolys,
                                                                slot, "polygons")), proj4string = CRS(proj4string(world.map))),
                                  world.map@data)

whole <- world.map[grep("Antarctica", world.map$NAME, invert=TRUE),]
half <- lSPDF[grep("Antarctica", lSPDF$NAME, invert=TRUE),]

whole <- fortify(whole, region="ISO3")
half <- fortify(half, region="ISO3")

world.map$scaled_growth <- as.numeric(scale(world.map@data$growth, 
                                            center = min(world.map@data$growth), 
                                            scale = max(world.map@data$growth)))

growth <- world.map@data[,c("ISO3", "scaled_growth")]
colnames(growth) <- c("id", "scaled_growth")
growth$scaled_growth <- factor(as.numeric(cut(growth$scaled_growth, 8))) # make it discrete

half_centers <- data.frame(cbind(coordinates(gCentroid(lSPDF, byid = TRUE)),
                                 id=world.map@data$ISO3, category=world.map@data$category))
half_centers$category <- factor(half_centers$category)

gg <- ggplot()
gg <- gg + geom_map(data=whole, map=whole, aes(x=long, y=lat, map_id=id), alpha=0, color="black", size=0.15)
gg <- gg + geom_map(data=growth, map=whole, aes(fill=scaled_growth, map_id=id))
gg <- gg + geom_map(data=half, map=half, aes(x=long, y=lat, map_id=id), fill="white")
gg <- gg + geom_point(data=half_centers, aes(x=x, y=y, shape=category), size=2)
gg <- gg + scale_fill_brewer(palette="Pastel2")
gg <- gg + scale_shape_discrete()
gg <- gg + coord_equal()
gg

enter image description here

like image 102
hrbrmstr Avatar answered Nov 05 '22 05:11

hrbrmstr


I think you can get (effectively) two different fill scales, with a little hack of scale_fill_brewer and scale_fill_manual.

Here's my output: enter image description here

I use the first bit of code from the other thread you posted in the question:

library(rgdal)
library(ggplot2)
library(maptools)

world.map <- readOGR(dsn="data", layer="TM_WORLD_BORDERS_SIMPL-0.3")

# Get centroids of countries
theCents <- coordinates(world.map)

# extract the polygons objects
pl <- slot(world.map, "polygons")

# Create square polygons that cover the east (left) half of each country's bbox
lpolys <- lapply(seq_along(pl), function(x) {
  lbox <- bbox(pl[[x]])
  lbox[1, 2] <- theCents[x, 1]
  Polygon(expand.grid(lbox[1,], lbox[2,])[c(1,3,4,2,1),])
})

# Slightly different data handling
wmRN <- row.names(world.map)

n <- nrow(world.map@data)
world.map@data[, c("growth", "category")] <- list(growth = 4*runif(n),
                                                  category = factor(sample(1:5, n, replace=TRUE)))

# Determine the intersection of each country with the respective "left polygon"
lPolys <- lapply(seq_along(lpolys), function(x) {
  curLPol <- SpatialPolygons(list(Polygons(lpolys[x], wmRN[x])),
                             proj4string=CRS(proj4string(world.map)))
  curPl <- SpatialPolygons(pl[x], proj4string=CRS(proj4string(world.map)))
  theInt <- gIntersection(curLPol, curPl, id = wmRN[x])
  theInt
})

# Create a SpatialPolygonDataFrame of the intersections
lSPDF <- SpatialPolygonsDataFrame(SpatialPolygons(
  unlist(lapply(lPolys,slot, "polygons")), 
  proj4string = CRS(proj4string(world.map))),
  world.map@data)

Now my contribution (borrowing the names whole/half from user hrbrmstr!)

# get two data.frames, one with whole countries and the other with the left half
# this relies on code from SO user BenBarnes
whole <- fortify(world.map, region="ISO3")
half <- fortify(lSPDF, region="ISO3")

# random growth / category data, similar to the random data originally
# suggested by Xu Wang
set.seed(123)
df <- data.frame(id = unique(world.map@data$ISO3),
                 growth = 4*runif(n),
                 category = factor(sample(letters[1:5], n, replace=T)))

# make growth a factor; 5 levels for convenience
df$growth_fac <- cut(df$growth, 5)

# append growth and category factor levels together
growth_cat_levels <- c(levels(df$category), levels(df$growth_fac))

# adjust factors with new joint levels
df$growth_fac <- 
  factor(df$growth_fac, levels=growth_cat_levels)
df$category <-
  factor(df$category, levels=growth_cat_levels)

# create a palette with some sequential colors and some qualitative colors
pal <- c(scale_fill_brewer(type='seq', palette=6)$palette(5),
         scale_fill_brewer(type='qual', palette='Pastel2')$palette(5))


# merge data
whole <- data.frame(merge(whole, df, by='id'))
half <- data.frame(merge(half, df, by='id'))

# plot
ggplot() +
  geom_polygon(data=whole, 
               aes(x=long, y=lat, group=group, fill=growth_fac), 
               color='black', size=0.15) +
  geom_polygon(data=half,
               aes(x=long, y=lat, group=group, fill=category), 
               color=NA) +
  scale_shape_discrete() +
  coord_equal() +

  scale_fill_manual('Category, Growth', 
                    values=pal, breaks=growth_cat_levels) +
  guides(fill=guide_legend(ncol=2))

Some notes:

  • I still think this is a hard-to-read map, but a fun challenge
  • I changed the 'category' names from numeric to alpha, to help avoid confusion with the 'growth' data.
  • I also kept the 'growth' data labels as generated by cut, to help make it clear that this is binned continuous data.
  • At first, I had the Growth colors on the left hand side of the legend, but I swapped it around; since Category determines the fill color of the left side country polygons, I thought Category should appear on the left in the legend
  • I experimented with a couple different palette options. One danger is that the qualitative scale will have a color too similar to the range of the sequential scale (as was the case with my post prior to this edit). Having one side greyscale and one side with colors helps avoid this
like image 23
arvi1000 Avatar answered Nov 05 '22 05:11

arvi1000