如何将 csv（栅格）读入 R 并将其与 shapefile 相交？

Question

我正在尝试将R中的csv数据集（csv_1）和

多边形形状文件

（cities_shp）放在一起。 csv 是分辨率为 100 米的栅格，其中有一个名为 “variable_1” 的感兴趣变量（如下图所示为多个点）。 shapefile 具有世界各地多个城市的多边形（如下图所示为单个多边形）。

我正在尝试创建 csv 文件和多边形 shapefile 的交集，以便我可以识别整个城市的“variable_1”的值。如下图所示，csv 点在地理上与“cities_shp”位置不对应。

#####################
# Inputs
#####################
## Libraries ##
library(terra)
library(sf)
library(dplyr)
library(mapview)
library(raster)

# Read City polygon
cities_shp = st_read("Input/.../GHS_STAT_UCDB2015MT_GLOBE_R2019A_V1_1.shp") 
cities_shp = st_make_valid(cities_shp)

# Read CSV file
csv_1 = read.csv("Input/.../csv_1.csv") 


#####################
# Graph 
#####################
csv_1_graph = st_as_sf(csv_1, coords = c("x", "y"), crs = 3395)
csv_1_graph = st_transform(csv_1_graph, crs = st_crs(cities_shp))
mapview(csv_1_graph) + mapview(cities_shp)

接下来我做了以下步骤：

1. 尝试光栅化“csv_1”文件，然后，
2. 将此新栅格转换为形状文件，
3. 使用步骤 2 中的 shapefile，这样我就可以将它与“cities_shp”相交。

我只能实现第二步，因为我无法使用适当的坐标来可视化这个新的形状文件。使用我使用的两个选项中的任何一个（如下），我不知道如何纠正 CRS，以便我可以将其与“cities_shp”相交。因为当我尝试时，下面显示的任何选项都会出现错误。

有谁知道如何读取这个 csv 文件，然后将其与 shapefile 相交？

注意：在本例中我只向您展示了一个城市，但我必须向全世界展示 10,000 个城市。因此，非常感谢任何帮助！

#####################
# 1st try
#####################
aux_2 <- try(rast(csv_1))
terra::plot(aux_2)
terra::plot(aux_2$variable_1)
aux_3 <- raster(aux_2)
mapview(aux_3,zcol = "variable_1",legend = T)
# Warning message:
#   In rasterCheckAdjustProjection(x, method) :
#   supplied layer has no projection information and is shown without background map
aux_2_polygons = as.polygons(aux_2)   |> st_as_sf()
class(aux_2_polygons)
new_df <- aux_2_polygons %>% st_as_sf(sf_column_name = "geometry", crs = 4326)
mapview(new_df)
#writing: substituting ENGCRS["Undefined Cartesian SRS with unknown unit"] for missing CRS

#####################
# 2nd try
#####################
geoName <- csv_1
geoName_raster <- raster(xmn = min(geoName$x), xmx = max(geoName$x), 
                         ymn = min(geoName$y), ymx = max(geoName$y), 
                         res = 100,crs = "+proj=longlat +datum=WGS84")
geoName_rasterize <- rasterize(geoName[, c('x', 'y')], geoName_raster, 
                               geoName[, 'variable_1'], fun = sum)
plot(geoName_rasterize)
points_matrix <- rasterToPoints(geoName_rasterize)
points_df <- as.data.frame(points_matrix) 
points_sf <- st_as_sf(points_df, coords = c("x", "y"), crs = "+proj=longlat +datum=WGS84")
points_sf = st_transform(points_sf, crs = 4326) #4326
plot(points_sf)
mapview(points_sf)

Answer 1

您为 csv_1 数据分配了错误的 CRS。 似乎正确的 CRS 是 WGS84 Pseudo Mercator/EPSG:3857。

但是，如果不知道 csv 的来源或它是如何创建的，这可能不正确。您需要请数据的创建者确认。

无论哪种方式：

library(sf)
library(mapview)

# Assign EPSG:3857 CRS
csv_1_graph <- st_as_sf(csv_1, coords = c("x", "y"), crs = 3857)

# Project to st_crs(cities_shp)
csv_1_graph <- st_transform(csv_1_graph, st_crs(cities_shp))

mapview(csv_1_graph) +
  mapview(cities_shp)