解决南美洲Shapefile与PM2.5栅格数据坐标不匹配叠加问题

我有包含南美洲Shapefile几何信息的数据框，尝试与同区域PM2.5栅格数据叠加，但两张地图无法对应，推测是数据排列或编码问题导致。

原代码 #

library(ncdf4)
library(lubridate)
library(RColorBrewer)
library(raster)
library(rasterVis)
library(chron)
library(ggplot2)
library(sf)

color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))(100)


shapefile_path <- "C:/Users/littl/South_America/South_America.shp"
south_america <- st_read(shapefile_path)
# 检查Shapefile的坐标系
st_crs(south_america)

# 如果Shapefile坐标系与地图不匹配，执行转换
# 例如，如果地图坐标系是WGS 84 (EPSG:4326)，可以执行：
if (st_crs(south_america)$epsg != 4326) {
south_america <- st_transform(south_america, crs = 4326)
}


# 确保输入正确的"yearmean_levtype_sfc.nc"文件路径
file_path <- "C:/Users/littl/yearmean_levtype_sfc.nc"
nc_file <- nc_open(file_path)
nc_file
# 提取第一波段数据(pm2p5)
pm2p5_data <- ncvar_get(nc_file, "pm2p5", start = c(1, 1, 1), count = c(-1, -1, 1)) 
#start = c(1, 1, 1) 指定从三个维度（经度、纬度、时间）的第一个位置开始读取数据。
#count = c(-1, -1, 1) 表示读取前两个维度（经度和纬度）的所有元素，仅读取时间维度的第一个波段。前两个维度的"-1"表示"读取所有元素"。
pm2p5_data
# 提取数据后关闭NetCDF文件
nc_close(nc_file)
# 转换为数据框
lon <- seq(-90, -20, length.out = dim(pm2p5_data)[1])
lat <- seq(20, -60, length.out = dim(pm2p5_data)[2])
df <- expand.grid(lon = lon, lat = lat)
#df$lat <- rev(df$lat)
#df$lon <- rev(df$lon)
df$pm2p5 <- as.vector(pm2p5_data)
ggplot() +
geom_raster(data = df, aes(x = lon, y = lat, fill = pm2p5)) +
scale_fill_gradientn(colors = color_palette) + # 使用自定义色板
geom_sf(data = south_america, fill = "transparent", color = "black", size = 0.5) +
labs(title = "Annual mean PM2.5, 2010",
     x = "经度",
     y = "纬度",
     fill = "PM2.5 (kg m^-3)") +
theme_minimal()

数据详情 #

South_America.shp信息 #

south_america
Simple feature collection with 15 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -109.4461 ymin: -58.49861 xmax: -26.24139 ymax: 12.59028
Geodetic CRS:  WGS 84
First 10 features:
                  COUNTRY                       geometry
1               Argentina MULTIPOLYGON (((-71.01648 -...
2                 Bolivia MULTIPOLYGON (((-62.19884 -...
3                  Brazil MULTIPOLYGON (((-70.62862 -...
4                   Chile MULTIPOLYGON (((-73.61806 -...
5                Colombia MULTIPOLYGON (((-81.71306 1...
6                 Ecuador MULTIPOLYGON (((-78.70903 -...
7  French Guiana (France) MULTIPOLYGON (((-54.60378 2...
8   Falkland Islands (UK) MULTIPOLYGON (((-59.79139 -...
9                  Guyana MULTIPOLYGON (((-60.08087 5...
10               Suriname MULTIPOLYGON (((-54.09445 3...

PM2.5 NetCDF文件信息 #

File C:/Users/littl/yearmean_levtype_sfc.nc (NC_FORMAT_64BIT):

     2 variables (excluding dimension variables):
        double time_bnds[bnds,time]   
        short pm2p5[longitude,latitude,time]   
            standard_name: mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air
            long_name: Particulate matter d <= 2.5 um
            units: kg m**-3
            add_offset: 2.38108791659124e-06
            scale_factor: 7.26694719096392e-11
            _FillValue: -32767
            missing_value: -32767
            cell_methods: time: mean

     4 dimensions:
        time  Size:10   *** is unlimited *** 
            standard_name: time
            long_name: time
            bounds: time_bnds
            units: hours since 1900-01-01 00:00:00.0
            calendar: gregorian
            axis: T
        bnds  Size:2 (no dimvar)
        longitude  Size:69 
            standard_name: longitude
            long_name: longitude
            units: degrees_east
            axis: X
        latitude  Size:94 
            standard_name: latitude
            long_name: latitude
            units: degrees_north
            axis: Y

    5 global attributes:
        CDI: Climate Data Interface version 2.0.4 (https://mpimet.mpg.de/cdi)
        Conventions: CF-1.6
        history: Wed Jul 19 14:23:12 2023: cdo yearmean levtype_sfc.nc yearmean_levtype_sfc.nc
2023-07-18 14:07:17 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/tmp/048068bd-b521-42fe-a5ff-0d5e1fe2fe60-adaptor.mars.internal-1689689093.9330037-19537-15-tmp.nc /cache/tmp/048068bd-b521-42fe-a5ff-0d5e1fe2fe60-adaptor.mars.internal-1689682802.5415428-19537-11-tmp.grib
        frequency: year
        CDO: Climate Data Operators version 2.0.4 (https://mpimet.mpg.de/cdo)

输出对比 #

• R中叠加效果：[错位的PM2.5栅格与南美洲边界]
• QGIS中叠加效果：[正确匹配的PM2.5栅格与南美洲边界]

问题根源并非坐标系不匹配（两者均为WGS84），而是栅格数据的纬度排列顺序与坐标映射不匹配，且未对NetCDF的编码数据进行解码。

修改步骤 #

1. 对NetCDF中的编码数据进行解码，还原真实PM2.5数值
2. 反转栅格数据的纬度维度，确保数据与坐标序列一一对应
3. 处理NetCDF中的填充值，避免异常值干扰绘图

修改后的代码 #

library(ncdf4)
library(lubridate)
library(RColorBrewer)
library(raster)
library(rasterVis)
library(chron)
library(ggplot2)
library(sf)

color_palette <- colorRampPalette(c("blue", "green", "yellow", "red"))(100)

# 读取Shapefile（已确认是WGS84，无需转换）
shapefile_path <- "C:/Users/littl/South_America/South_America.shp"
south_america <- st_read(shapefile_path)

# 读取并处理NetCDF数据
file_path <- "C:/Users/littl/yearmean_levtype_sfc.nc"
nc_file <- nc_open(file_path)

# 读取编码后的PM2.5数据及元数据
pm2p5_encoded <- ncvar_get(nc_file, "pm2p5", start = c(1,1,1), count = c(-1,-1,1))
scale_factor <- ncatt_get(nc_file, "pm2p5", "scale_factor")$value
add_offset <- ncatt_get(nc_file, "pm2p5", "add_offset")$value
fill_value <- ncatt_get(nc_file, "pm2p5", "_FillValue")$value

# 解码数据并替换填充值为NA
pm2p5_data <- pm2p5_encoded * scale_factor + add_offset
pm2p5_data[pm2p5_encoded == fill_value] <- NA

# 反转纬度维度，匹配坐标序列
pm2p5_data <- pm2p5_data[, ncol(pm2p5_data):1]

nc_close(nc_file)

# 创建坐标数据框
lon <- seq(-90, -20, length.out = dim(pm2p5_data)[1])
lat <- seq(20, -60, length.out = dim(pm2p5_data)[2])
df <- expand.grid(lon = lon, lat = lat)
df$pm2p5 <- as.vector(pm2p5_data)

# 绘图
ggplot() +
  geom_raster(data = df, aes(x = lon, y = lat, fill = pm2p5)) +
  scale_fill_gradientn(colors = color_palette, na.value = "transparent") +
  geom_sf(data = south_america, fill = "transparent", color = "black", size = 0.5) +
  labs(title = "2010年PM2.5年均值",
       x = "经度",
       y = "纬度",
       fill = "PM2.5 (kg m^-3)") +
  theme_minimal()

关键说明 #

• 数据解码：NetCDF中的pm2p5是经过缩放编码的短整型，必须用scale_factor和add_offset计算得到真实值，否则数值完全错误
• 纬度反转：NetCDF数据按从北到南存储纬度，但as.vector按列展开数据，反转纬度维度后才能让数据与坐标一一对应
• 填充值处理：将NetCDF中的无效填充值替换为NA，避免绘图时出现异常色块

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