как отрезать или обрезать или заполнить белым цветом большой. расширенный (на 10%) прямоугольник за пределами многоугольника с ggplot2

этот вопрос является продолжением моего предыдущий вопрос SO и связан с этот вопрос.

я просто пытаюсь заполнить белым цветом область на 10% больше, чем простой многоугольник с помощью ggplot2. может я неправильно группирую? вот фото шипа с воспроизводимым кодом ниже

введите здесь описание изображения

# 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

r gis ggplot2 map map-projections
person Anthony Damico    schedule 14.10.2014    source источник
comment
Из документации функции coord_map: Это все еще экспериментально, и если у вас есть какие-либо советы относительно лучшего (или более правильного) способа сделать это, пожалуйста, дайте мне знать. Могу поспорить, что вы нашли ошибку в этой экспериментальной функции.   -  person Pop    schedule 17.10.2014
comment
Вы имеете в виду на 10% больше, чем размеры ограничивающей рамки полигона?   -  person jbaums    schedule 17.10.2014
comment
@jbaums да, знаю. так что для этого примера 10% за пределами самой дальней границы штата Коннектикут во всех направлениях :)   -  person Anthony Damico    schedule 17.10.2014
comment
Вы настроены на ggplot2 решение? Предполагая, что я понимаю ваши намерения, этот материал довольно прост с базовым графиком.   -  person jbaums    schedule 17.10.2014
comment
Чтобы упростить пример кода, вы можете заменить все строки, участвующие в буферизации и различении ваших фигур, строками 5-7 из этого недавнего вопроса. . Я также должен согласиться с jbaums: вчера я запустил ваш код и проголосовал за вопрос, но как только я понял, что это проблема с ggplot и coord_map в частности, я перестал копать, так как ggplot чертовски сложно (для меня, по крайней мере) отлаживать! С функцией построения базы или решетки я бы, по крайней мере, попробовал.   -  person Josh O'Brien    schedule 17.10.2014
comment
@jbaums, я надеюсь решить эту проблему с помощью ggplot2, потому что другие компоненты полного кода, который я пытаюсь усовершенствовать полагаться на geom_tile(). тем не менее, я, безусловно, был бы признателен за решение, отличное от ggplot2, если вы считаете его простым   -  person Anthony Damico    schedule 17.10.2014
comment
@JoshO'Brien очень подлый! я отредактировал свой пример и буду использовать его в будущем. спасибо!   -  person Anthony Damico    schedule 17.10.2014
comment
@ ДжошО'Брайен и Энтони: осторожно - 1.1 * extent(x) - это не то же самое, что добавить 10% к каждой стороне, если только обе оси не охватывают экватор. В данном случае 1.1 * xmax и 1.1 * ymax фактически вычитают 10% из этих границ (поскольку xmax и ymax отрицательны).   -  person jbaums    schedule 17.10.2014
comment
@jbaums -- Не совсем так. raster определяет метод для * применительно к объекту класса Extent, и он умнее этого. Попробуйте extent(10,20,10,20) * 1.2, чтобы убедиться, что он действительно дополняет каждое измерение. (Обратите внимание, что вам нужно умножить на 1.2, чтобы получить 10% расширение в каждом направлении.) Если вы хотите взглянуть на код, лежащий в основе метода * Extent, вы можете получить его, набрав getMethod("Arith", c("Extent", "numeric")).   -  person Josh O'Brien    schedule 17.10.2014
comment
@JoshO'Brien, ах, спасибо за объяснение - я не знал о таком поведении! Действительно, очень хитро :)   -  person jbaums    schedule 17.10.2014

Ответы (2)


arrow_upward
12
arrow_downward

Это связано с тем, что `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)

введите здесь описание изображения

person kohske    schedule 19.10.2014
comment
(+1) Хорошо. Не могли бы вы объяснить, что делают две строки с hoge? Я не вижу, чтобы is.linear.hoge использовалось снова позже (или это вызывается coord_map?). - person jbaums; 19.10.2014
comment
Вот это да. как и @jbaums, мне интересно узнать больше о том, что здесь происходит и зачем нужны эти дополнительные строки ... но это явно решает представленную проблему. благодарю вас!! - person Anthony Damico; 19.10.2014
comment
немного поиграл с этим, похоже, что включение coord_fixed() в ваш последний вызов необходимо, если ваши слои включают что-то, что не является geom_polygon - person Anthony Damico; 19.10.2014

arrow_upward
5
arrow_downward

Вот как я поступил бы с базовыми функциями построения графиков. Мне было не совсем ясно, нужно ли, чтобы «фоновый» полигон отличался от полигона состояния, или можно ли, чтобы он был простым прямоугольником, на который будет наложен полигон состояния. Возможно и то, и другое, но я сделаю последнее здесь для краткости/простоты.

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)

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

[EDIT: я внес некоторые изменения в следующий код. Теперь он (необязательно) отображает линии сетки, что особенно важно при построении осей в единицах, которые находятся в другой проекции графика.]

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)

карта

person jbaums    schedule 17.10.2014