在 R 中处理不规则间隔的网格 netcdf 数据答案

【问题标题】：Working with irregularly spaced gridded netcdf data in R在 R 中处理不规则间隔的网格 netcdf 数据
【发布时间】：2021-05-08 20:23:54
【问题描述】：

我很快就会处理模拟的空气污染数据。我得到了一些示例数据，看起来网格非常不规则（'brick' 会抛出关于不规则网格的错误），这使得将其转换为栅格格式变得更加困难。

    # load packages ####
    library(ncdf4)
    library(raster)
    library(rgdal)
    library(fields)
    
    # netcdf to data matrix ####
    nc_data <- nc_open('filename.nc')

    cnc_PM2_5       <- ncvar_get(nc_data, attributes(nc_data$var)$names[1])
    longitude       <- ncvar_get(nc_data, "longitude")
    latitude        <- ncvar_get(nc_data, "latitude")
    #time            <- ncvar_get(nc_data, "time")

    # plot matrix, latitude is upside down ####
    image.plot(longitude,rev(latitude), PM2_5, 
               main="PM2_5", ylab="latitude")

    # get quick shapefile for Finland (http://www.diva-gis.org/datadown)####
    finland <- readOGR("FIN_adm4.shp")
    # add map projection
    proj4string(finland) <- CRS("+proj=longlat +datum=WGS84 +init=EPSG:4326")
    helsinki <- finland[finland$NAME_4=="Helsinki",]
    # add Helsinki outline to image.plot
    plot(helsinki, add=TRUE)

如何将这些不规则数据放入空间对象中？（SpatialPointsDataFrame 或 raster）

可重现的例子

    PM25 <- matrix(data=c(4.540674,3.437939,3.373072,2.335204,2.140603,2.529804,
                          4.475807,4.346074,6.42181,3.113605,2.594671,2.854138,
                          3.632539,6.097476,3.308205,5.643409,5.643409,2.140603,
                          4.151473,5.643409,2.529804,2.400071,3.827139,4.994741,
                          3.243339,3.373072,2.075737,2.400071,4.281207,3.697406,
                          2.854138,5.448809,2.075737,3.827139,3.567672,2.854138,
                          4.475807,2.983871,3.697406,5.578542,2.140603,2.20547,
                          6.22721,3.502806,5.773142,4.994741,5.254208,2.594671,
                          3.502806,3.827139,3.502806,2.400071,3.762273,2.724404,
                          3.892006,3.502806,2.724404,2.20547,4.151473,5.838009,
                          6.746144,5.059608,2.140603,2.270337,3.567672,2.140603,
                          3.437939,2.724404,2.140603,3.437939,3.308205,4.994741,
                          4.086606,2.335204,2.140603,5.578542,4.994741,2.724404,
                          3.892006,4.605541,4.281207,4.281207,2.140603,2.529804,
                          3.308205,3.892006,3.697406,3.308205,2.270337,6.097476,
                          3.308205,4.41094,5.319075,4.086606,2.854138,6.162343,
                          2.140603,3.308205,3.048738,4.02174), ncol=10, nrow=10)
    
    # latitude
    latitude <- c(60.1191177368164,60.1192321777344,60.1193504333496,
                          60.1194686889648,60.1195831298828,60.119701385498,
                          60.1198196411133,60.1199340820312,60.1200523376465,
                          60.1201667785645)
    # longitude
    longitude <- c(24.5949993133545,24.595235824585,24.5954704284668,
                          24.5957050323486,24.5959415435791,24.5961761474609,
                          24.5964126586914,24.5966472625732,24.5968818664551,
                          24.5971183776855)
    
    image.plot(longitude,latitude, PM25, main="PM2_5")

【问题讨论】：

标签： r geospatial netcdf

【解决方案1】：

只是一个快速更新，使用stars 包和read_ncdf 更容易。

library(stars)
AQI <- read_ncdf('filename.nc', var = c("cnc_PM2_5")

【讨论】：

【解决方案2】：

感谢@Spacedman 提供this 答案：

将存储为矩阵的数据转换为坐标点的向量：

先把坐标展开成一整套坐标点：

latlong = expand.grid(long=longitude, lat=latitude)

然后通过使用 c(.) 展平矩阵来添加包含数据的列：

latlong = cbind(latlong, PM25 = c(PM25))

给予：

head(latlong)
      long      lat PM25
1 24.59500 60.11912    1
2 24.59524 60.11912    2
3 24.59547 60.11912    3
4 24.59571 60.11912    4
5 24.59594 60.11912    5
6 24.59618 60.11912    6

然后用sp设置坐标：

coordinates(latlong)=~long+lat

绘制时，它应该看起来与 image.plot 相同，除了点和可能不同的颜色：

spplot(latlong,"PM25")

【讨论】：