wie man schnippt oder beschneidet oder ein großes weiß füllt. erweitert (um 10%) Rechteck außerhalb eines Polygons mit ggplot2

Question

Diese Frage ist ein Follow-up meiner prior SO question und ist verwandt mit this question.

Ich versuche nur einen Bereich zu füllen, der 10% größer ist als ein einfaches Polygon mit ggplot2. Vielleicht gruppiere ich Dinge falsch? hier ist ein Foto von der Spitze mit reproduzierbaren Code unter

# 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 

Answer 1

Dies liegt daran, `coord_map‘, oder allgemeinere nichtlineare Koordinaten interpoliert intern Scheitelpunkte, so dass die Linie als eine der Koordinate entsprechende Kurve gezeichnet wird. In Ihrem Fall wird die Interpolation zwischen einem Punkt des äußeren Rechtecks ​​und einem Punkt der inneren Kante durchgeführt, den Sie als Bruch sehen.

Sie können dies ändern, indem Sie:

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

Sie können auch den Unterschied des Verhaltens finden Sie hier:

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) 

Answer 2

Hier ist, wie ich es mit Basis Plotfunktionen gehen würde. Es war mir nicht ganz klar, ob das "Hintergrund" -Polygon Unterschiede zu dem Zustandspolygon benötigt oder ob es ein einfaches Rechteck ist, bei dem der Zustand poly überlagert ist. Entweder ist möglich, aber ich mache das letztere hier aus Gründen der Kürze/Einfachheit.

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) 

Nun, wenn ich richtig verstehe, obwohl sie in einem projizierten Koordinatensystems ist, wollen Sie Achsen zeichnen, die in den Einheiten eines anderen sind (das Original) Koordinatensystem. Hier ist eine Funktion, die das für Sie tun kann.

[BEARBEITEN: Ich habe einige Änderungen am folgenden Code vorgenommen. Es ist nun (optional) zeichnet die Gitterlinien, die besonders wichtig sind, wenn Achse in Einheiten darstellen, die in einem anderen Vorsprung an dem Grundstück sind.]

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)