使用 ggplot 在 R 中绘制组合 shapefile 时出现问题，未给出定义错误

Posted 2023-02-16

【中文标题】使用 ggplot 在 R 中绘制组合 shapefile 时出现问题，未给出定义错误

【英文标题】：Issues plotting a combined shapefile in R using ggplot, no defined error given

【发布时间】：2020-10-23 07:55:17

【问题描述】：

使用的库：

library(sp)
library(sf) 
library(ggplot2) 
library(ggmap)

创建了一个名为“coordinate.data”的数据框，其中经度和纬度作为列名，气象站位置作为行名。

longitude <- c(-73.964482,-73.953678,-73.893522,-73.815856,-74.148499)
latitude <- c(40.767544,40.631762,40.872481,40.734335,40.604014)
coordinate.data <- data.frame(longitude,latitude)
rownames(coordinate.data) <- c("MANH","BKLN","BRON","QUEE","STAT")

然后我检索了新泽西县和纽约市行政区的 shapefile 数据，并删除了所有不必要的列，因此两个 shapefile 中只剩下几何字段。 NYC Boroughs shapefile 数据是从 NYC Open Data 下载的，而 NJ 县边界是从 NJGIN Open Data 下载的。

nj.shp <- st_read("~/Downloads/NJ/NJ_Counties.shp")
nj <- nj.shp[,-(1:21)]
nyc.shp <- st_read("~/Downloads/NY/NYC_Boroughs.shp")
nyc <- nyc.shp[,-(1:4)]

我将两个 shapefile 格式化为具有相同的投影（ESPG 代码 3857）并将它们组合成一个 shapefile 数据框，其中包含一个变量（几何）中的 26 个观测值（县/区）。

same.projection <- CRS("+init=EPSG:3857")
nj.data <- st_transform(nj,same.projection)
new.projection <- CRS("+init=EPSG:3857")
nyc.data <- st_transform(nyc,new.projection)
combined.data <- rbind(nj.data,nyc.data)

我现在正尝试在地图上绘制组合 shapefile（“combined.data”）以及气象站位置（“coordinate.data”）。 当我尝试这样做时，它不可避免地会运行并且 R 会关闭。如果我删除 geom_sf(...)，它会正确绘制站点并格式化所有内容，所以我认为问题出在这行代码上。

mesonet.map <-ggplot() +
  ggtitle("NY Mesonet Site Locations") +
  xlab("Longitude") +
  ylab("Latitude") +
  geom_point(data=coordinate.data,aes(x=longitude,y=latitude))+
  geom_text(aes(x=longitude,y=latitude,label=rownames(coordinate.data)),size=3.25,nudge_y=0.02)+
  geom_sf(data=combined.data,fill='darkgreen') +

mesonet.map + theme(
  panel.background=element_rect(fill="lightblue",color="lightblue",size=0.5,linetype="solid"),
  panel.grid.major=element_line(size=0.5,linetype='solid',color="white"), 
  panel.grid.minor=element_line(size=0.25,linetype='solid',color="white")
)

【问题讨论】：

除了一些小问题，您的代码在我的机器上运行良好。

【参考方案1】：

我不太确定您从这些形状文件中删除了什么。我没有删除任何东西。我也没有结合任何东西。在我的输出中投影了四个单独的图层：基本地图、NY 形状文件、NJ 形状文件和称为站点的形状文件，该文件在下面创建和描述。

此地图是使用 3 个现有形状文件和在以下步骤中创建的站点形状文件创建的。从每个步骤创建的变量都会打印出来，并与每个步骤的描述一起显示。

首先，创建空间几何变量：

MULTIPOINT ((-73.96448 40.76754), (-73.95368 40.63176), (-73.89352
 40.87248), (-73.81586 40.73434), (-74.1485 40.60401))

然后，创建几何列（结合多点变量和crs）

MULTIPOINT ((-73.96448 40.76754), (-73.95368 40...

Geometry set for 1 feature 
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: -74.1485 ymin: 40.60401 xmax: -73.81586 ymax: 40.87248
projected CRS:  WGS 84 / Pseudo-Mercator

然后创建一个由站点名称或一列组成的常规 df。

然后创建了简单的特征对象 - （将几何列与 df 相结合）。输出显示：

A sf: 5 × 2
site    pts.sfc
<fct>   <MULTIPOINT [m]>
MANH    MULTIPOINT ((-73.96448 40.7...
BKLN    MULTIPOINT ((-73.96448 40.7...
BRON    MULTIPOINT ((-73.96448 40.7...
QUEE    MULTIPOINT ((-73.96448 40.7...
STAT    MULTIPOINT ((-73.96448 40.7...

然后这个 sf 对象被写入 shape 文件夹：

Writing layer `Weather_sites' to data source `C:/pathR' using driver `ESRI Shapefile'
Writing 5 features with 1 fields and geometry type Multi Point.

输出显示纽约州以及新泽西州和纽约州的 5 个地点。 NY 形状文件似乎不那么具有包容性，仅包含一些有限的 NY 数据。

正确完成整个映射需要按顺序执行四个步骤。几何，几何列，sf对象，写入shape文件夹。

The map exported to .png