2012-11-04 73 views
20

我使用的例子在這裏討論: 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") 

得出以下結果: enter image description here

我想映射一個變量一個國家的左半部分,另一個國家的右半部分。我把「half」放在引號中,因爲它沒有明確定義(或者至少我沒有明確定義它)。 Ian Fellows的回答可能會有所幫助(這可以讓您輕鬆獲得質心)。我希望的東西,這樣我可以在例如做aes(left_half_color = growth, right_half_color = category)。如果情況不同,我也對上半部和下半部感興趣。

如果可能的話,我也想了半個人重心映射到一些東西。

+7

您可能要考慮有兩個並排的地圖。可能比這個國家的分裂更直觀地看待和解釋。 –

+0

@Marcinthebox謝謝你的建議。 –

回答

26

這是一個沒有ggplot一個解決方案,依賴於plot函數來代替。它還要求rgeos包除了在OP的代碼:

EDIT現在有10%以下的視覺疼痛

編輯2現在提供質心爲東西半部

library(rgeos) 
library(RColorBrewer) 

# 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([email protected]) 
[email protected][, 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))), 
    [email protected]) 

########## 
## EDIT ## 
########## 
# Create a slightly less harsh color set 
s_growth <- scale([email protected]$growth, 
    center = min([email protected]$growth), scale = max([email protected]$growth)) 
growthRGB <- colorRamp(c("red", "blue"))(s_growth) 
growthCols <- apply(growthRGB, 1, function(x) rgb(x[1], x[2], x[3], 
    maxColorValue = 255)) 
catCols <- brewer.pal(nlevels([email protected]$category), "Pastel2") 

# and plot 
plot(world.map, col = growthCols, bg = "grey90") 

plot(lSPDF, col = catCols[[email protected]$category], add = TRUE) 

enter image description here

也許有人能想出w^ith一個很好的解決方案,使用ggplot2。然而,基於this answer到有關多個填充尺度單個圖形問題(「你不能」),一個ggplot2解決方案似乎不大可能無刻面(這可能是一個很好的方法,因爲在意見提出以上)。


編輯回覆:半的東西映射重心:的質心爲東(「左」)可以得到一半的

coordinates(lSPDF) 

者爲西(「右」 )半部可以通過以類似的方式創建rSPDF對象獲得:

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

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

# Create a SpatialPolygonDataFrame of the western (right) intersections 
rSPDF <- SpatialPolygonsDataFrame(SpatialPolygons(unlist(lapply(rPolys, 
    slot, "polygons")), proj4string = CRS(proj4string(world.map))), 
    [email protected]) 

然後信息可以在地圖上根據繪製lSPDFrSPDF的質心:

points(coordinates(rSPDF), col = factor([email protected]$REGION)) 
# or 
text(coordinates(lSPDF), labels = [email protected]$FIPS, cex = .7) 
+0

謝謝你的好回答(和更新)。如果我按照以下網站上的建議,將允許我結合自己做過的事,但對於GGPLOT2? https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles –

+0

@XuWang,你應該能夠使用鏈接的指令來繪製'lSPDF'和'rSPDF' shape文件,但據我所知,你會遇到的問題如果你想爲每個半部分填充不同的映射。 – BenBarnes

+0

感謝您的幫助和回覆/更新。 –