географические карты тепла с использованием решетки levelplot в R

Question

Я пытаюсь создать геопространственную непрерывную карту тепла, которая не ограничена округами или искусственными границами, используя решетку, очень похожую на Dialect Maps от Katz, и немного отличается от 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), вероятно, хорошо работать, но я не вижу простой способ сделать это. Есть предположения?

Answer 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") 

Answer 2

Вы можете преобразовать 'pts' и «m» к объектам класса Raster* и Spatial* соответственно, а затем использовать их как вход для spplot (который является оберткой вокруг levelplot, см. ?raster::spplot). Это позволяет отбрасывать точно те пиксели, которые не должны отображаться (через mask), и в то же время обеспечивает совместимость с дальнейшими решеткой (или). Вот пример кода.

## 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.