spTransfrom 函数的替代方案由于 rgdal 的退役

Question

我有一个脚本可以使用 shapefile 从 gbif 中过滤出现数据。由于 rgdal/rgeos 退役，我必须做出一些改变。

head(sp_gbif)
        lon       lat     gbifID
1 -70.95618 -53.27092 3457120752
2 -68.31865 -54.81420 3031584233
3  22.76239  66.60015 3431642207
4  16.69931  64.20202 3432283190
5  23.32176  66.11418 3432398188
6  22.29740  66.74461 3432398189

shape
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 7.511393, 10.49182, 47.5338, 49.79137  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 6
names       : USE, RS,       RS_ALT,                GEN,       SHAPE_LENG,       SHAPE_AREA 
value       :   2, 08, 080000000000, Baden-Württemberg, 1298891.66320858, 35801397076.4735


#subset the GBIF data into a data frame
occ.map <- data.frame(sp_gbif[["lon"]], sp_gbif[["lat"]], sp_gbif[["gbifID"]])
print(str(occ.map, 1))

#simplify column names
names(occ.map)[1:3] <- c('long', 'lat', 'gbifID')
print(head(occ.map, 10))

#turning the data frame into a "spatial points data frame"
coordinates(occ.map) <- c("long", "lat")

#defining the datum 
proj4string(occ.map) <- crs("+proj=longlat +datum=WGS84")
crs(shape) <- "+proj=longlat +datum=WGS84"

#reprojecting the 'gbif' data frame to the same as in the 'shape' object

occ.map <- spTransform(occ.map, proj4string(shape))

#Identifying records from gbif that fall within the shape polygons
inside <- occ.map[apply(gIntersects(occ.map, shape, byid = TRUE), 2, any),]

#Prepare data frame for joining with the occcurrence df so only records 
#that fall inside the polygons get selected 
res.gbif <- data.frame(inside@data)
final.gbif <- sp_gbif %>% semi_join(res.gbif, by = "gbifID")
head(final.gbif)

我尝试使用 sf 但 sf_transform 例如不适用于空间点数据帧。

有人有调整代码的想法吗？

Answer 1

为了使用

st_transform

包的

sf

函数将 lon-lat 坐标转换为笛卡尔（平面）坐标，您必须切换到使用

sf

空间数据框。首先，将您的

sp

空间数据框转换为具有以下功能的

sf

空间数据框：

vz_sf <- st_as_sf(vz)    #to convert from `sp` spatial data frame or;
vz_sf <- st_read(vz.shp) #to read from standard data formats like Shapefiles

注意：如果可以的话，要读取文件，请使用 GeoPackage。

一个大问题可能是，如果您最终需要使用

sp

数据框调用包中的函数，以及使用

sf

数据框的包中的函数。在这种情况下，您应该尝试尽可能多地使用

sf

来保持工作流程，然后在需要互操作性时从一种格式转换为另一种格式：

sfdata = st_as_sf(spdata)
spdata = as(sfdata, "Spatial")