无法读取 netcdf 文件并将其转换为 SpatRaster

Trouble reading a netcdf file and convert it into a SpatRaster

我正在努力打开 NetCDF 文件并使用 R 将其转换 为栅格。 数据应该位于 25 公里乘 25 公里的规则网格上。它包含海 北极的冰浓度。

library(terra)
#> terra 1.5.21
library(ncdf4)

file <- "~/Downloads/data_sat_Phil_changt_grid/SIC_SMMR_month_2015.nc"

我收到有关未找到范围的警告。

r <- rast(file)
#> Error in R_nc4_open: No such file or directory
#> Warning: [rast] GDAL did not find an extent. Cells not equally spaced?

我们可以看到coordinates/extent有问题。

r
#> class       : SpatRaster 
#> dimensions  : 448, 304, 12  (nrow, ncol, nlyr)
#> resolution  : 0.003289474, 0.002232143  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : SIC_SMMR_month_2015.nc:sic 
#> varname     : sic 
#> names       : sic_1, sic_2, sic_3, sic_4, sic_5, sic_6, ...

我可以用 nc_open() 打开 nc 文件,我看到坐标存在。

nc <- nc_open(file)
names(nc$var)
#> [1] "lat" "lon" "sic"

lat <- ncvar_get(nc, "lat")
lon <- ncvar_get(nc, "lon")

dim(lat)
#> [1] 304 448
dim(lon)
#> [1] 304 448
dim(r)
#> [1] 448 304  12

是否可以 assemble 此数据(SIC 值和坐标)创建 SpatRaster? nc文件可以在这里下载:https://easyupload.io/pfth0s

reprex package (v2.0.1)

于 2022-05-20 创建

数据已网格化,但文件未指定坐标,也未指定坐标参考系。该文件指定了与单元格关联的 lon/lat 值,但对我们帮助不大,因为这些值显然不在规则网格上。从 plot(r)

很容易看出
NAflag(r) = -9999
plot(r,1)

也来自

p = cbind(as.vector(lon), as.vector(lat))
plot(p, cex=.1, xlab="lon", ylab="lat")

因此,您需要找出使用的是哪个坐标参考系 (crs),显然是某种极地 crs。以及数据集的范围是多少。

从你指向的网站来看,我认为我们可以使用:

crs(r) = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"
ext(r) = c(-3850000, 3750000, -5350000, 5850000)
r
#class       : SpatRaster 
#dimensions  : 448, 304, 12  (nrow, ncol, nlyr)
#resolution  : 25000, 25000  (x, y)
#extent      : -3850000, 3750000, -5350000, 5850000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs 
#source      : SIC_SMMR_month_2015.nc:sic 
#varname     : sic 
#names       : sic_1, sic_2, sic_3, sic_4, sic_5, sic_6, ... 
 

结果看起来不错:

g = geodata::gadm("Greenland", level=0, path=".")
gg = project(g, crs(r))
plot(r,1)
lines(gg)

但这当然不是做这种事情的好方法; ncdf 文件应该包含所有需要的元数据。