在R中读取CSV格式栅格数据并与Shapefile进行空间交集分析的方法

嘿，我看你现在的核心问题是CSV格式的栅格点数据和城市多边形Shapefile空间不匹配，还需要完成两者的空间交集来提取城市范围内的variable_1值对吧？别担心，我结合你已经提到的terra和sf库，一步步帮你搞定：

一、先解决空间不匹配的核心问题：统一坐标系 #

你说CSV点和Shapefile地理位置不对应，大概率是两个数据的坐标系不一样，这可是空间分析里最常见的坑！先从这里入手：

1. 把CSV数据转换成带空间参考的点对象 #

首先得把你的CSV数据转换成sf的空间点对象，前提是你的CSV里有经度、纬度或者其他空间坐标字段（比如x、y）：

# 读取CSV文件（替换成你的实际文件名）
csv_raw <- read.csv("your_csv_file.csv")
# 转换为sf空间点对象，coords参数填你的坐标字段名，crs填数据实际的坐标系（比如WGS84是4326）
csv_sf <- st_as_sf(csv_raw, coords = c("lon", "lat"), crs = 4326)

2. 读取Shapefile并确认坐标系 #

接着读取你的城市Shapefile，同时查看它的坐标系：

# 读取城市Shapefile（替换成你的实际文件名）
cities_shp <- st_read("your_cities_shapefile.shp")
# 查看Shapefile的坐标系，记下来这个CRS代码
st_crs(cities_shp)

3. 统一两个数据的坐标系 #

把CSV的点对象转换成和Shapefile完全一致的坐标系，这样空间位置就能对应上了：

# 转换CSV点的坐标系，和Shapefile保持同步
csv_sf_aligned <- st_transform(csv_sf, crs = st_crs(cities_shp))

二、执行空间交集分析 #

现在两个数据的空间基准统一了，就可以提取城市范围内的variable_1值了，分两种常用场景：

场景1：保留所有落在城市内的点数据，关联城市属性 #

如果你想拿到每个落在城市里的CSV点，同时带上对应的城市信息，用st_intersection就可以：

# 执行空间交集，结果会包含城市的所有属性和CSV里的variable_1
city_point_intersection <- st_intersection(csv_sf_aligned, cities_shp)
# 查看结果
head(city_point_intersection)

场景2：按城市多边形汇总variable_1的值 #

如果需要对每个城市的variable_1做统计（比如平均值、总和），可以先做空间连接再分组计算：

# 先把CSV点和城市多边形做空间连接，给每个点匹配所属的城市
city_variable_join <- st_join(csv_sf_aligned, cities_shp)
# 按城市分组统计，这里的city_name替换成你Shapefile里的城市名称字段
city_variable_summary <- city_variable_join %>%
  dplyr::group_by(city_name) %>%
  dplyr::summarise(
    avg_variable1 = mean(variable_1, na.rm = TRUE),
    total_variable1 = sum(variable_1, na.rm = TRUE),
    count_points = n()
  )

三、如果你的CSV是真正的栅格数据（不是点） #

你提到CSV是100米分辨率的栅格，要是它其实是ASCII栅格格式（比如包含行列号、分辨率等栅格元数据），用terra处理会更高效：

# 用terra读取栅格数据
raster_data <- rast("your_csv_raster.csv")
# 把sf的Shapefile转换成terra的矢量对象
cities_vect <- vect(cities_shp)
# 提取每个城市范围内的栅格值，还可以直接指定统计函数
city_raster_stats <- extract(raster_data, cities_vect, fun = mean, na.rm = TRUE)

备注：内容来源于stack exchange，提问作者Adriana Castillo Castillo