Карта Choropleth в ggplot с полигонами с отверстиями

Я пытаюсь нарисовать карту хороплета Германии, показывающую уровень бедности по штату (вдохновленный этим вопросом).

Проблема в том, что некоторые из государств (например, Берлин) полностью окружены другими государствами (Бранденбург), и мне трудно получить ggplot, чтобы узнать "дыру" в Бранденбурге.

Данные для этого примера здесь.

library(rgdal)
library(ggplot2)
library(RColorBrewer)

map <- readOGR(dsn=".", layer="germany3")
pov <- read.csv("gerpoverty.csv")

mrg.df <- data.frame(id=rownames([email protected]),[email protected]$ID_1)
mrg.df <- merge(mrg.df,pov, by="ID_1")
map.df <- fortify(map)
map.df <- merge(map.df,mrg.df[,c("id","poverty")], by="id")
ggplot(map.df, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(fill=poverty))+
  geom_path(colour="grey50")+
  scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

Обратите внимание, что цвета для Берлина и Бранденбурга (на северо-востоке) идентичны. Их не должно быть - уровень бедности в Берлине намного ниже, чем у Бранденбурга. Похоже, что ggplot отображает берлинский многоугольник, а затем передает над ним Бранденбургский полигон без отверстия.

Если я изменю вызов geom_polygon(...), как предложено здесь, я могу исправить проблему Берлина/Бранденбурга, но теперь три самых северных государства отображаются неправильно.

ggplot(map.df, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(group=poverty, fill=poverty))+
  geom_path(colour="grey50")+
  scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

Что я делаю неправильно?

Ответ 1

Вы можете построить многоугольники острова в отдельном слое, следуя примеру ggplot2 wiki. Я изменил шаги слияния, чтобы сделать это проще:

mrg.df <- data.frame(id=rownames([email protected]),[email protected]$ID_1)
mrg.df <- merge(mrg.df,pov, by="ID_1")
map.df <- fortify(map)
map.df <- merge(map.df,mrg.df, by="id")

ggplot(map.df, aes(x=long, y=lat, group=group)) +
    geom_polygon(aes(fill=poverty), color = "grey50", data =subset(map.df, !Id1 %in% c("Berlin", "Bremen")))+
    geom_polygon(aes(fill=poverty), color = "grey50", data =subset(map.df, Id1 %in%  c("Berlin", "Bremen")))+
    scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
    labs(x="",y="")+ theme_bw()+
    coord_fixed()

Как незапрашиваемый акт евангелизации, я рекомендую вам рассмотреть что-то вроде

library(ggmap)
qmap("germany", zoom = 6) +
    geom_polygon(aes(x=long, y=lat, group=group, fill=poverty),
                 color = "grey50", alpha = .7,
                 data =subset(map.df, !Id1 %in% c("Berlin", "Bremen")))+
    geom_polygon(aes(x=long, y=lat, group=group, fill=poverty),
                 color = "grey50", alpha= .7,
                 data =subset(map.df, Id1 %in%  c("Berlin", "Bremen")))+
    scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))

чтобы обеспечить контекст и знакомые ориентиры.

Ответ 2

Это просто расширение на ответ @Ista, которое не требует, чтобы кто-нибудь знал, какие состояния (Берлин, Бремен) нужно сделать последним.

В этом подходе используется тот факт, что fortify(...) генерирует столбец hole, который идентифицирует, является ли группа координат дырой. Таким образом, это отображает все области (id) с любыми отверстиями до (например, внизу) областей без отверстий.

Большое спасибо @Ista, без ответа я не мог бы придумать это (поверьте, я потратил много часов, пытаясь...)

ggplot(map.df, aes(x=long, y=lat, group=group)) +
  geom_polygon(data=map.df[map.df$id %in% map.df[map.df$hole,]$id,],aes(fill=poverty))+
  geom_polygon(data=map.df[!map.df$id %in% map.df[map.df$hole,]$id,],aes(fill=poverty))+
  geom_path(colour="grey50")+
  scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

Ответ 3

Просто добавьте еще одно небольшое улучшение в ответы @Ista и @jhoward (большое спасибо за вашу помощь!).

Модификация @jhoward может быть легко завернута в небольшую функцию, подобную этой

gghole <- function(fort){
        poly <- fort[fort$id %in% fort[fort$hole,]$id,]
        hole <- fort[!fort$id %in% fort[fort$hole,]$id,]
        out <- list(poly,hole)
        names(out) <- c('poly','hole')
        return(out)
} 
# input has to be a fortified data.frame

Тогда не нужно каждый раз вспоминать, как извлекать информацию о отверстиях. Код будет выглядеть как

    ggplot(map.df, aes(x=long, y=lat, group=group)) +
            geom_polygon(data=gghole(map.df)[[1]],aes(fill=poverty),colour="grey50")+
            geom_polygon(data=gghole(map.df)[[2]],aes(fill=poverty),colour="grey50")+
    # (optionally). Call by name
    #         geom_polygon(data=gghole(map.df)$poly,aes(fill=poverty),colour="grey50")+
    #         geom_polygon(data=gghole(map.df)$hole,aes(fill=poverty),colour="grey50")+
            scale_fill_gradientn(colours=brewer.pal(5,"OrRd"))+
            labs(x="",y="")+ theme_bw()+
            coord_fixed()

Ответ 4

В качестве альтернативы вы можете создать эту карту, используя rworldmap.

library(rworldmap)
library(RColorBrewer)
library(rgdal)

map <- readOGR(dsn=".", layer="germany3")
pov <- read.csv("gerpoverty.csv")

#join data to the map
sPDF <- joinData2Map(pov,nameMap='map',nameJoinIDMap='VARNAME_1',nameJoinColumnData='Id1')

#default map
#mapPolys(sPDF,nameColumnToPlot='poverty')

colours=brewer.pal(5,"OrRd")
mapParams <- mapPolys( sPDF
                      ,nameColumnToPlot='poverty'
                      ,catMethod="pretty"
                      ,numCats=5
                      ,colourPalette=colours
                      ,addLegend=FALSE )


do.call( addMapLegend, c( mapParams
                          , legendLabels="all"
                          , legendWidth=0.5
                        ))

#to test state names
#text(pov$x,pov$y,labels=pov$Id1)