栅格与矢量数据CRS转换及缓冲区提取技术咨询

数据CRS检查结果 #

> st_crs(cru_precip)
    Coordinate Reference System:
      User input: +proj=longlat +datum=WGS84 
      wkt:
    GEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]],
        CS[ellipsoidal,2],
            AXIS["longitude",east,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            AXIS["latitude",north,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]]]
> st_crs(districts)
    Coordinate Reference System:
      User input: WGS 84 
      wkt:
    GEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        CS[ellipsoidal,2],
            AXIS["latitude",north,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433]],
            AXIS["longitude",east,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]]

问题分析 #

两个数据均基于WGS84椭球的大地坐标系，核心差异是：

• cru_precip 对应EPSG:6326，轴序为经度在前、纬度在后
• districts 对应EPSG:4326，轴序为纬度在前、经度在后

这种轴序差异会导致空间匹配错位，且地理坐标系的度单位不适合创建实际距离的缓冲区（比如米/千米），必须针对性处理。

解决方案步骤 #

1. 统一地理坐标系，修正轴序
   优先将cru_precip转换为通用的EPSG:4326，确保两者空间对齐：

   cru_precip <- st_transform(cru_precip, 4326)
   

2. 转换到投影坐标系创建缓冲区
   地理坐标系的度单位无法准确表示实际距离，必须切换到投影坐标系（如UTM分带、区域专属投影）来生成缓冲区：

   • 选择适配研究区域的投影（例如中国区域可选用CGCS2000 UTM分带，如EPSG:4523；全球小区域用对应UTM分带）
   • 转换矢量数据、计算质心与缓冲区：

     # 示例：假设研究区位于UTM 50N（EPSG:32650），请根据实际区域调整
     districts_proj <- st_transform(districts, 32650)
     districts_centroids <- st_centroid(districts_proj)
     # 创建5000米缓冲区，按需修改dist参数值
     districts_buffers <- st_buffer(districts_centroids, dist = 5000)
     

3. 提取栅格数据
   两种可选方式，根据数据规模选择：

   • 方式一：缓冲区转回WGS84提取原栅格（适合小数据量）

     districts_buffers_wgs84 <- st_transform(districts_buffers, 4326)
     precip_extracted <- raster::extract(cru_precip, districts_buffers_wgs84)
     

   • 方式二：栅格转换到投影坐标系提取（精度更一致，适合大数据量）

     cru_precip_proj <- raster::projectRaster(cru_precip, crs = st_crs(districts_proj)$wkt)
     precip_extracted <- raster::extract(cru_precip_proj, districts_buffers)
     

关键提示 #

• 不要直接在地理坐标系下创建缓冲区，否则得到的是"度"单位的范围，不符合实际距离需求
• 投影坐标系的选择必须贴合研究区域，避免投影变形影响缓冲区精度

内容的提问来源于stack exchange，提问作者crawling_ec