如何剪断或裁剪或填充大的。使用 ggplot2 在多边形外扩展（增加 10%）矩形答案

【问题标题】：how to snip or crop or white-fill a large. expanded (by 10%) rectangle outside of a polygon with ggplot2如何剪断或裁剪或填充大的。使用 ggplot2 在多边形外扩展（增加 10%）矩形
【发布时间】：2014-10-14 11:38:24
【问题描述】：

这个问题是我prior SO question的后续问题，与this question有关。

我只是想用 ggplot2 填充比简单多边形大 10% 的区域。也许我把事情分组错了？这是尖峰的照片，下面有可重现的代码

# reproducible example
library(rgeos)
library(maptools)
library(raster)

shpct.tf <- tempfile() ; td <- tempdir()

download.file( 
    "ftp://ftp2.census.gov/geo/pvs/tiger2010st/09_Connecticut/09/tl_2010_09_state10.zip" ,
    shpct.tf ,
    mode = 'wb'
)

shpct.uz <- unzip( shpct.tf , exdir = td )

# read in connecticut
ct.shp <- readShapePoly( shpct.uz[ grep( 'shp$' , shpct.uz ) ] )

# box outside of connecticut
ct.shp.env <- gEnvelope( ct.shp )
ct.shp.out <- as( 1.2 * extent( ct.shp ), "SpatialPolygons" )


# difference between connecticut and its box
ct.shp.env.diff <- gDifference( ct.shp.env , ct.shp )
ct.shp.out.diff <- gDifference( ct.shp.out , ct.shp )


library(ggplot2)


# prepare both shapes for ggplot2
f.ct.shp <- fortify( ct.shp )
env <- fortify( ct.shp.env.diff )
outside <- fortify( ct.shp.out.diff )


# create all layers + projections
plot <- ggplot(data = f.ct.shp, aes(x = long, y = lat))  #start with the base-plot 

layer1 <- geom_polygon(data=f.ct.shp, aes(x=long,y=lat), fill='black')

layer2 <- geom_polygon(data=env, aes(x=long,y=lat,group=group), fill='white')

layer3 <- geom_polygon(data=outside, aes(x=long,y=lat,group=id), fill='white')

co <- coord_map( project = "albers" , lat0 = 40.9836 , lat1 = 42.05014 )

# this works
plot + layer1 

# this works
plot + layer2

# this works
plot + layer1 + layer2

# this works
plot + layer2 + co

# this works
plot + layer1 + layer3

# here's the problem: this breaks
plot + layer3 + co

# this also breaks, but it's ultimately how i want to display things
plot + layer1 + layer3 + co

# this looks okay in this example but
# does not work for what i'm trying to do-
# cover up points outside of the state
plot + layer3 + layer1 + co

【问题讨论】：

来自coord_map 函数的文档：“这仍然是实验性的，如果您有任何关于更好（或更正确）方法的建议，请告诉我” .我敢打赌，你在这个“实验”功能中发现了一个错误
你的意思是比多边形的边界框的尺寸大10%吗？
@jbaums 是的，我愿意。所以对于这个例子，在各个方向上都超出了康涅狄格州边界的最远范围 10% :)
为了简化您的示例代码，您可以用来自this recent question 的第 5-7 行替换所有涉及缓冲和区分形状的行。我也必须同意 jbaums：我昨天运行了你的代码并支持了这个问题，但是一旦我意识到这是 ggplot 和coord_map 的问题，我就停止了挖掘，因为 ggplot 太难了（至少对我来说) 调试！使用 base 或 lattice 绘图功能，我至少可以试一试。
@jbaums -- 不太正确。 raster 为* 定义了一个方法，当应用于Extent 类对象时，它比这更聪明。试试extent(10,20,10,20) * 1.2 看看它确实添加到每个维度。（请注意，您必须乘以 1.2 才能在每个方向上获得 10% 的扩展。）如果您想查看 Extent 的 * 方法的底层代码，您可以通过键入getMethod("Arith", c("Extent", "numeric")).

标签： r map ggplot2 gis map-projections

【解决方案1】：

这是因为`coord_map'，或更一般的非线性坐标，在内部对顶点进行插值，以便将线绘制为对应于坐标的曲线。在您的情况下，将在外部矩形的点和内部边缘的点之间执行插值，您将其视为中断。

您可以通过以下方式更改：

co2 <- co
class(co2) <- c("hoge", class(co2))
is.linear.hoge <- function(coord) TRUE
plot + layer1 + layer3 + co2

你也可以在这里找到行为的区别：

ggplot(data.frame(x = c(0, 90), y = 45), aes(x, y)) + geom_line() + co + ylim(0, 90)

ggplot(data.frame(x = c(0, 90), y = 45), aes(x, y)) + geom_line() + co2 + ylim(0, 90)

【讨论】：

(+1) 不错。你介意解释一下hoge 的两行代码在做什么吗？我没有看到is.linear.hoge 以后再次使用（或者它是由coord_map 调用的？）。
哇。像@jbaums 一样，我很好奇这里发生了什么以及为什么这些额外的行是必要的。但这显然解决了所提出的问题。谢谢！！
玩了一下，如果你的层包含任何不是geom_polygon的东西，看起来在你的最终调用中包含coord_fixed()是必要的

【解决方案2】：

以下是我将如何使用基本绘图功能进行处理。我并不完全清楚你是否需要“背景”多边形与状态多边形不同，或者它是否可以是一个简单的矩形，将覆盖状态多边形。任何一种都是可能的，但为了简洁起见，我会在这里做后者。

library(rgdal)
library(raster) # for extent() and crs() convenience

# download, unzip, and read in shapefile
download.file(file.path('ftp://ftp2.census.gov/geo/pvs/tiger2010st/09_Connecticut/09',
                        'tl_2010_09_state10.zip'), f <- tempfile(), mode='wb')
unzip(f, exdir=tempdir())
ct <- readOGR(tempdir(), 'tl_2010_09_state10')

# define albers and project ct
# I've set the standard parallels inwards from the latitudinal limits by one sixth of
# the latitudinal range, and the central meridian to the mid-longitude. Lat of origin 
# is arbitrary since we transform it back to longlat anyway.
alb <- CRS('+proj=aea +lat_1=41.13422 +lat_2=41.86731 +lat_0=0 +lon_0=-72.75751
           +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')
ct.albers <- spTransform(ct, alb)

# expand bbox by 10% and make a polygon of this extent
buf <- as(1.2 * extent(ct.albers), 'SpatialPolygons')
proj4string(buf) <- alb

# plot without axes
par(mar=c(6, 5, 1, 1)) # space for axis labels
plot(buf, col='white', border=NA)
do.call(rect, as.list(c(par('usr')[c(1, 3, 2, 4)], col='gray90'))) 
# the above line is just in case you needed the grey bg
plot(buf, add=TRUE, col='white', border=NA) # add the buffer
plot(ct.albers, add=TRUE, col='gray90', border=NA)
title(xlab='Longitude')
title(ylab='Latitude', line=4)

现在，如果我理解正确，尽管处于投影坐标系中，但您想要绘制以另一个（原始）坐标系为单位的轴。这是一个可以为您执行此操作的函数。

[编辑：我对以下代码进行了一些更改。它现在（可选）绘制网格线，这在以与绘图不同投影的单位绘制轴时尤其重要。]

axis.crs <- function(plotCRS, axisCRS, grid=TRUE, lty=1, col='gray', ...) {
  require(sp)
  require(raster)
  e <- as(extent(par('usr')), 'SpatialPolygons')
  proj4string(e) <- plotCRS
  e.ax <- spTransform(e, axisCRS)
  if(isTRUE(grid)) lines(spTransform(gridlines(e.ax), plotCRS), lty=lty, col=col)
  axis(1, coordinates(spTransform(gridat(e.ax), plotCRS))[gridat(e.ax)$pos==1, 1],
       parse(text=gridat(e.ax)$labels[gridat(e.ax)$pos==1]), ...)
  axis(2, coordinates(spTransform(gridat(e.ax), plotCRS))[gridat(e.ax)$pos==2, 2],
       parse(text=gridat(e.ax)$labels[gridat(e.ax)$pos==2]), las=1, ...)
  box(lend=2) # to deal with cases where axes have been plotted over the original box
}

axis.crs(alb, crs(ct), cex.axis=0.8, lty=3)

【讨论】：