GGPLOT2/GIS多邊形區域

Question

內繪製的再生用，請考慮以下數據：

library(rgdal) 
library(ggplot2) 
library(rgeos) 

download.file("http://spatialanalysis.co.uk/wp-content/uploads/2010/09/London_Sport.zip", 

destfile = "London_Sport.zip") 
unzip("London_Sport.zip") 

projection="+proj=merc" 

london_shape = readOGR("./", layer="london_sport") 

# Create random points 
set.seed(1) 
points = data.frame(long=rnorm(10000, mean=-0.1, sd=0.1), lat=rnorm(10000, mean=51.5, sd=0.1)) 
points = SpatialPoints(points, proj4string=CRS("+proj=latlon")) 

# Transform data to our projection 
london = spTransform(london_shape, CRS(projection)) 
points = spTransform(points, CRS(projection)) 

# Keeps only points inside London 
intersection = gIntersects(points, london, byid = T) 
outside = apply(intersection == FALSE, MARGIN = 2, all) 
points = points[which(!outside), ] 

# Blank theme 
new_theme_empty <- theme_bw() 
new_theme_empty$line <- element_blank() 
new_theme_empty$rect <- element_blank() 
new_theme_empty$strip.text <- element_blank() 
new_theme_empty$axis.text <- element_blank() 
new_theme_empty$plot.title <- element_blank() 
new_theme_empty$axis.title <- element_blank() 
new_theme_empty$plot.margin <- structure(c(0, 0, -1, -1), unit = "lines", valid.unit = 3L, class = "unit") 

# Prepare data to ggplot 
london = fortify(london) 
points = as.data.frame(points) 

我想積點的密度圖。我可以通過使用stat_bin2d這樣做：

ggplot() + 
    geom_polygon(data=london, aes(x=long,y=lat,group=group), fill="black") + 
    stat_bin2d(data=points, aes(x=long,y=lat), bins=40) + 
    geom_path(data=london, aes(x=long,y=lat,group=id), colour='white') + 
    coord_equal() + 
    new_theme_empty 

但是，這導致密度廣場的某些部分可以被繪製倫敦之外：

我怎麼能只繪製密度圖在倫敦？

Answer 1

我找到了答案，得到的是倫敦的邊界框和倫敦本身（使用gDifference）之間的區別的多邊形，並將其繪製成白色之上。我認爲這種方法的缺點是：1）如果某些方塊仍然出現在它後面，則必須手動增加多邊形的尺寸。 2）你不能使用複雜背景的繪圖主題。因此，如果任何人有更好的答案，我會留下一段時間的問題。

下面的代碼：

library(rgdal) 
library(ggplot2) 
library(rgeos) 

projection="+proj=merc" 

#London boroughs polygons 
download.file("http://spatialanalysis.co.uk/wp-content/uploads/2010/09/London_Sport.zip", destfile = "London_Sport.zip") 
unzip("London_Sport.zip") 
london = readOGR("./", layer="london_sport") 
london = spTransform(london, CRS(projection)) 

# Generate random points 
set.seed(1) 
points = data.frame(long=rnorm(10000, mean=-0.1, sd=0.1), lat=rnorm(10000, mean=51.5, sd=0.1)) 
points = SpatialPoints(points, proj4string=CRS("+proj=latlon")) 
points = spTransform(points, CRS(projection)) 

# Keep only points inside London 
intersection = gIntersects(points, london, byid = TRUE) 
inside = apply(intersection == TRUE, MARGIN = 2, any) 
points = points[which(inside), ] 

# Create a bounding box 10% bigger than the bounding box of London 
x_excess = (london@bbox['x','max'] - london@bbox['x','min'])*0.1 
y_excess = (london@bbox['y','max'] - london@bbox['y','min'])*0.1 
x_min = london@bbox['x','min'] - x_excess 
x_max = london@bbox['x','max'] + x_excess 
y_min = london@bbox['y','min'] - y_excess 
y_max = london@bbox['y','max'] + y_excess 
bbox = matrix(c(x_min,x_max,x_max,x_min,x_min, 
       y_min,y_min,y_max,y_max,y_min), 
       nrow = 5, ncol =2) 
bbox = Polygon(bbox, hole=FALSE) 
bbox = Polygons(list(bbox), "bbox") 
bbox = SpatialPolygons(Srl=list(bbox), pO=1:1, proj4string=london@proj4string) 

# Get the Polygon that is the difference between the bounding box and London 
outside = gDifference(bbox,london) 

# Blank theme 
new_theme_empty <- theme_bw() 
new_theme_empty$line <- element_blank() 
new_theme_empty$rect <- element_blank() 
new_theme_empty$strip.text <- element_blank() 
new_theme_empty$axis.text <- element_blank() 
new_theme_empty$plot.title <- element_blank() 
new_theme_empty$axis.title <- element_blank() 
new_theme_empty$plot.margin <- structure(c(0, 0, -1, -1), unit = "lines", valid.unit = 3L, class = "unit") 

# Prepare data for ggplot 
london = fortify(london) 
points = as.data.frame(points) 
outside = fortify(outside) 

# Plot! 
ggplot() + 
    geom_polygon(data=london, aes(x=long,y=lat,group=group), fill="black") + 
    stat_bin2d(data=points, aes(x=long,y=lat), bins=40) + 
    geom_path(data=london, aes(x=long,y=lat,group=id), colour='white') + 
    geom_polygon(data=outside, aes(x=long,y=lat), fill='white') + 
    coord_equal() + 
    new_theme_empty