在 R 中使用晶格水平图的地理热图答案

【问题标题】：geographical heat maps using lattice levelplot in R在 R 中使用晶格水平图的地理热图
【发布时间】：2014-09-24 14:56:28
【问题描述】：

我正在尝试使用格子构建一个不受县或人工边界限制的地理空间连续热图，很像 Katz 的 Dialect Maps，与 Choropleth Challenge 有点不同。

到目前为止，我已经很接近了，但我一直在努力弄清楚如何避免在地理区域以外的区域显示热图颜色。请参见下面的示例：

library(maps)
library(lattice)

# get region border coordinates for the contiguous USA
m <- map("usa")

# make a grid of latitude and longitude, and supply z-values
lons <- seq(min(m$x, na.rm=T), max(m$x, na.rm=T), length.out=30)
lats <- seq(min(m$y, na.rm=T), max(m$y, na.rm=T), length.out=30)
pts <- expand.grid(lons, lats)
names(pts) <- c("lon", "lat")
pts$z <- sin(pts$lat*pi/180) + cos(pts$lon*pi/180)

## (A) eliminate z-values outside of the USA region
# pts$z[ lat & lon outside of region ] <- NA # don't know how to do this
# this could work, but would leave jagged edges around the region border

levelplot(z~lon*lat, pts, aspect="xy",
           panel = function(...){
             panel.levelplot(...) 
             panel.xyplot(m$x, m$y, type="l", col="black") # adds USA border
             # (B) fill the area outside region with white 
             # panel.something() # not sure what to use here
          })

方法 (B) 可能会很好，但我没有看到一个简单的方法来做到这一点。有什么想法吗？

【问题讨论】：

这可能会有所帮助：r-sig-geo.2731867.n2.nabble.com/…

标签： r map heatmap lattice

【解决方案1】：

我不确定如何使用map 和lattice 完成您的要求，但我将使用raster 和rgeos 工具解决您的问题：

library(raster)
library(rgeos)

## get SpatialPolygnsDataFrame map of the states
m <- getData("GADM", country="United States", level=1)
m <- m[!m$NAME_1 %in% c("Alaska","Hawaii"),] # sorry Alaska and Hawaii 

## here I modified your code to make a raster object
r <- raster(nrow=30, ncol=30, 
            xmn=bbox(m)["x","min"], xmx=bbox(m)["x","max"],
            ymn=bbox(m)["y","min"], ymx=bbox(m)["y","max"],
            crs=proj4string(m))
xyz <- rasterToPoints(r)
r[] <- sin(xyz[,"y"]*pi/180) + cos(xyz[,"x"]*pi/180)

## Option A) mask raster using polygon
newr <- mask(r, m)
plot(newr, col=cm.colors(60), axes=FALSE)
plot(m, add=TRUE)
box(col="white")
## leaves jagged edges...

## Option B) cover the outside area
b <- gUnaryUnion(rasterToPolygons(r)) # first create a polygon that covers the raster
b <- gDifference(b, m) # then get the difference between the polygons
plot(r, col=cm.colors(100), axes=FALSE)
plot(b, add=TRUE, col="white", border="white")
plot(m, add=TRUE)
box(col="white")

【讨论】：

【解决方案2】：

您可以将 'pts' 和 'm' 分别转换为 Raster* 和 Spatial* 类的对象，然后将它们用作 spplot 的输入（这是 levelplot 的包装器，请参阅 @987654327 @)。这允许您准确丢弃那些不应该显示的像素（通过mask），同时确保与基于lattice（或grid）的进一步兼容操作。这是一些示例代码。

## transform 'map' object to 'SpatialPolygons'
library(maptools)
m <- map2SpatialPolygons(m, IDs = seq(m$names), 
                         proj4string = CRS("+init=epsg:4326"))

## rasterize pts and mask pixels outside map region
library(raster)
coordinates(pts) <- ~ lon + lat
proj4string(pts) <- "+init=epsg:4326"
pts <- as(pts, "SpatialPixelsDataFrame")
rst <- raster(pts)
rst <- mask(rst, m)

## display data
library(latticeExtra)
spplot(rst, scales = list(draw = TRUE), alpha.regions = .8) + 
  layer(sp.polygons(m, lwd = 2))

为了适当解决边缘效应，我建议您坚持使用@PaulRegular 提供的不错的解决方法。

【讨论】：