繪圖顏色編碼的網格和多邊形的交集

Question

我有一個多邊形（科羅拉多的輪廓）和覆蓋多邊形的顏色編碼的緯度 - 經度網格。我希望僅在多邊形本身內顯示顏色編碼網格。我怎樣才能做到這一點？

繪製多邊形和網格的代碼如下。爲了說明，我已將多邊形放置在下圖中的網格上方，否則您將無法看到多邊形。

我試過使用PBSmapping包來獲得兩個多邊形的交集。但是，我只能使用單一顏色繪製該交叉區域，而不是讓一個多邊形的原始顏色方案保留。 （網格是多色的。）

我也想過嘗試讓多邊形的內部沒有填充，而多邊形以外的區域是白色的。也許這是可能的。然後我可以將多邊形放在網格上，網格只能在多邊形內可見。但是，如果我有多個多邊形和一個網格，那麼這種方法可能不起作用。

謝謝你的任何建議。這是迄今爲止代碼：

library(rgdal) 
library(maptools) 
library(RColorBrewer) 
library(classInt) 
library(raster) 

setwd('c:/users/mark w miller/gis_in_R') 

## Specify a geographic extent for the map 
## by defining the top-left and bottom-right geographic coordinates 
mapExtent <- rbind(c(-115, 43), c(-97, 35)) 

# assign projection 
newProj <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") 

## Project the map extent (first need to specify that it is longlat) 
mapExtentPr <- spTransform(SpatialPoints(mapExtent, 
        proj4string=CRS("+proj=longlat")), 
        newProj) 

# get Colorado polygon 
us1 <- getData('GADM', country="USA", level=1) 
colorado <- us1[(us1$NAME_1 %in% c('Colorado')),] 

## Project Colorada layer 
coloradoPr <- spTransform(colorado, newProj) 

# create grid of polygons 
grd <- GridTopology(c(-110.5, 35.5), c(1,1), c(11,8)) 
polys <- as(grd, "SpatialPolygons") 

# assign projection to grid 
proj4string(polys) = CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") 

# create fake atttribute data for each grid cell 
poly.data = data.frame(f=runif(length(row.names(polys)), 0, 14)) 
row.names(poly.data) <- paste0('g', 1:length(row.names(polys))) 

# convert grid to a SpatialPolygonsDataFrame: 
poly.df = SpatialPolygonsDataFrame(polys, poly.data) 

# assign colors to grid cells 
plotvar <- poly.df$f 
nclr <- 8 
plotclr <- brewer.pal(nclr,"BuPu") 

colcode <- ifelse((    plotvar <= 2), plotclr[1], 
      ifelse((plotvar > 2 & plotvar <= 4), plotclr[2], 
      ifelse((plotvar > 4 & plotvar <= 6), plotclr[3], 
      ifelse((plotvar > 6 & plotvar <= 8), plotclr[4], 
      ifelse((plotvar > 8 & plotvar <= 10), plotclr[5], 
      ifelse((plotvar > 10 & plotvar <= 12), plotclr[6], 
      ifelse((plotvar > 12 & plotvar <= 14), plotclr[7], 
                plotclr[8]))))))) 

jpeg(filename = "colorado.with.grid2.jpeg") 

## Plot each projected layer, beginning with the projected extent 
plot(mapExtentPr, pch=NA) 
plot(poly.df, , col=colcode, add=TRUE) 
plot(coloradoPr , border="white", col="lightgrey", add=TRUE) 

dev.off() 

編輯

一個可能不精確，很費力的解決方案可能是創建多個網格是結合時大約每填多邊形。這與黎曼和會有點類似。但是，我必須認爲有一個更準確和易於實施的解決方案。

Answer 1

以下是我想出了到目前爲止，使用PBSmapping包：

library(rgdal) 
library(maptools) 
library(RColorBrewer) 
library(classInt) 
library(raster) 

setwd('c:/users/mark w miller/gis_in_R') 

## Specify a geographic extent for the map 
my.map.extent <- matrix(c(-113, 46, 
          -113, 33, 
          -93, 33, 
          -93, 46, 
          -113, 46), nrow=5, byrow=TRUE) 

my.map.extent.p = Polygon(my.map.extent) 
my.map.extent.ps = Polygons(list(my.map.extent.p),1) 
my.map.extent.sps = SpatialPolygons(list(my.map.extent.ps)) 
proj4string(my.map.extent.sps) <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") 

# assign projection 
newProj <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") 

# get my.states polygon 
us1 <- getData('GADM', country="USA", level=1) 
my.states <- us1[(us1$NAME_1 %in% c('Wyoming', 'Kansas')),] 

## Project state layer 
my.statesPr <- spTransform(my.states, newProj) 

# create grid of polygons 
#      lower left, cell size, number of cells on x axis and on y axis 
grd <- GridTopology(c(-112.5, 33.5), c(1,1), c(20,15)) 
polys <- as(grd, "SpatialPolygons") 

# assign projection to grid 
proj4string(polys) = CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") 

# create fake attribute data for each grid cell 
poly.data = data.frame(f=runif(length(row.names(polys)), 0, 14)) 
row.names(poly.data) <- paste0('g', 1:length(row.names(polys))) 

# convert grid to a SpatialPolygonsDataFrame: 
poly.df = SpatialPolygonsDataFrame(polys, poly.data) 

# assign colors to grid cells 
plotvar <- poly.df$f 
nclr <- 8 
plotclr <- brewer.pal(nclr,"BuPu") 

colcode <- ifelse((    plotvar <= 2), plotclr[1], 
      ifelse((plotvar > 2 & plotvar <= 4), plotclr[2], 
      ifelse((plotvar > 4 & plotvar <= 6), plotclr[3], 
      ifelse((plotvar > 6 & plotvar <= 8), plotclr[4], 
      ifelse((plotvar > 8 & plotvar <= 10), plotclr[5], 
      ifelse((plotvar > 10 & plotvar <= 12), plotclr[6], 
      ifelse((plotvar > 12 & plotvar <= 14), plotclr[7], 
                plotclr[8]))))))) 

library(PBSmapping) 

extent.ps <- combinePolys(SpatialPolygons2PolySet(my.map.extent.sps)) 
grid.ps <- combinePolys(SpatialPolygons2PolySet(poly.df)) 
states.ps <- combinePolys(SpatialPolygons2PolySet(my.statesPr)) 

diff.grid.states <- combinePolys(joinPolys(extent.ps, states.ps, "DIFF")) 

dev.new() 

jpeg(filename = "my.states.with.grid.jpeg") 

plotPolys(extent.ps, proj = TRUE, xlab="Easting", ylab="Northing") 
plot(poly.df, col=colcode, lwd = 0.5, add=TRUE) 
addPolys(diff.grid.states, col="white") 

dev.off()