如何使用levelplot或spplot生成状态shapefile和空间点？

Question

我有一个光栅图像，我想要与夏威夷的shapefile叠加，以便为地图提供更好的上下文。此外，我想在完成的地图上绘制特定的点。目前，我有一个夏威夷shapefile的ggplot图像，其中绘制了点：

df <- data.frame(lon, lat)

ggplot(data = world) + geom_sf() + 
       coord_sf(xlim = c(-160, -150), ylim = c(10, 25)) +
       geom_point(data = df, mapping = aes(x = lon, y = lat) +
       theme()

我的栅格是使用levelplot构建的

library(raster)
b <- brick(system.file("external/rlogo.grd", package = "raster"))
levelplot(b)

当我尝试：

> ggplot(data = world) + geom_sf() + coord_sf(xlim = c(-158, -154), ylim = c(17, 24), expand = FALSE) + levelplot(b)

我得到了这样的信息：

Error in ggplot(data = world) + geom_sf() + coord_sf(xlim = c(-158, -154),  : 
  non-numeric argument to binary operator
In addition: Warning message:
Incompatible methods ("+.gg", "+.trellis") for "+" 

如果这些映射类型不兼容，我如何通过levelplot或spplot生成包含空间点的夏威夷shapefile？我意识到通过使用可复制的Rstudio光栅，这些数据是古怪的，但我想让这个问题变得可重现。

提前感谢您的指导！

Answer 1

你可以从this网站下载夏威夷的shapefile。然后，您可以使用以下代码绘制它

library(sf)
library(ggplot2)
library(rgdal)
library(rgeos)

#Reading the shapefiles
sf <- st_read(dsn="C:\\Users\\User\\Desktop\\cty_council_dist.shp", layer="cty_council_dist_haw")
shape <- readOGR(dsn="C:\\Users\\User\\Desktop\\cty_council_dist.shp", layer="cty_council_dist_haw")

#To view the attributes
head(shape@data)
summary(sf)

#Plotting the shapefile
plot(shape)
plot(sf)

#Covert the coordinate system
sf_gcs <- st_transform(sf, crs = "EPSG:4326")
shape_gcs <- spTransform(shape, CRS=CRS("+init=EPSG:4326"))

#Plotting the shapefile
plot(shape)
plot(sf_gcs)

#Plotting the districts only
plot(sf_gcs["cntydist"], axes = TRUE, main = "Districts")

​

现在添加您可以使用的点数据

df <- data.frame(lon = c(-155.5, -155.2, -155.3), lat = c(19.2, 19.4, 19.8))

由于您尚未提供点数据，因此我创建了一些假数据

#使用ggplot2绘图

ggplot() + 
 geom_sf(data = sf_gcs, aes(fill = cntydist)) + theme(legend.position = "none") +
 geom_point(df, mapping = aes(x = lon, y = lat), col="green")

​

更新

要使用ggplot2绘制栅格，可以使用以下代码

#Load the libraries
library(rasterVis)
library(RColorBrewer)

## Create a matrix with random data & use image()
xy <- matrix(rnorm(400),20,20)
image(xy)

# Turn the matrix into a raster
rast <- raster(xy)

# Give it lat/lon coords for 156.2-154.8°W, 18.5-20.5°N
extent(rast) <- c(-156.2, -154.8, 18.5, 20.5)

#Assign a projection
projection(rast) <- CRS("+proj=longlat +datum=WGS84")
plot(rast)

colr <- colorRampPalette(brewer.pal(11, 'RdYlBu'))

#Plotting Using ggplot2
gplot(rast) +  
  geom_tile(aes(fill=factor(value),alpha=0.8)) + 
  geom_polygon(data=shape_gcs, aes(x=long, y=lat, group=group),
               fill=NA,color="grey50", size=1) + 
  theme(legend.position = "none")

​

​

如果您想使用rasterVis包中的levelplot绘制它，您可以使用

levelplot(rast, 
          margin=F,                       # suppress marginal graphics
          colorkey=list(
            space='right'                   # plot legend at right
          ),    
          par.settings=list(
            axis.line=list(col='transparent') # suppress axes and legend outline
          ),
          scales=list(draw=FALSE),            # suppress axis labels
          col.regions=colr,                   # colour ramp
          at=seq(-5, 5, len=101)) +           # colour ramp breaks
  layer(sp.polygons(shape_gcs, lwd=1))           # add shapefile with latticeExtra::layer

​

​