从 netcdf 到栅格 (tif) 的转换会改变分辨率和范围？答案

【问题标题】：Conversion from netcdf to raster (tif) changes resolution and extent?从 netcdf 到栅格 (tif) 的转换会改变分辨率和范围？
【发布时间】：2022-08-03 06:02:10
【问题描述】：

我正在尝试将 netcdf 文件转换为光栅 (tif) 格式。我创建了一个脚本，不久前它运行良好。但是现在，当我尝试对不同的文件使用相同的简单脚本时，分辨率从0.5 x 0.5 变为0.5 x 0.5263158。范围也从：

-100.25, -73.25, 28.75, 48.75

至

-100.5, -73, 28.48684, 49.01316

我也尝试过在 R 中使用不同的光栅包，但它们返回一条消息说单元格不等间距。文件（附加here）可能是一个很好的问题，但我看不出在哪里以及如何。

复制代码：


# load netcdf file
import xarray as xr 
import rioxarray

xds = xr.open_dataset(\'output_shocks_us/hybrid_gfdl-esm4_ssp126_2015co2_yield_soybean_shift_2017-2044.nc\')
xds = xds.rename({\'lat\':\'y\',\'lon\':\'x\', \'time\':\'band\'})

# Add CRS
xds.rio.write_crs(\"epsg:4326\", inplace=True)

# Convert to geotiff
xds[\"yield-soy-noirr\"].rio.to_raster(\'hybrid_gfdl-esm4_ssp126_2015co2_yield_soybean_shift_2017-2044_test.tif\')
rio = xr.open_rasterio(\"hybrid_gfdl-esm4_ssp126_2015co2_yield_soybean_shift_2017-2044_test.tif\")

print(xds)
print(rio)

完整的结果是：

print(xds)
<xarray.Dataset>
Dimensions:          (y: 39, x: 55)
Coordinates:
    band             int64 2025
  * y                (y) float64 28.75 30.25 30.75 31.25 ... 47.75 48.25 48.75
  * x                (x) float64 -100.2 -99.75 -99.25 ... -74.25 -73.75 -73.25
    spatial_ref      int32 0
Data variables:
    yield-soy-noirr  (y, x) float64 nan nan nan nan nan ... nan nan nan nan nan
Attributes:
    grid_mapping:  spatial_ref

############
print(rio)
<xarray.DataArray (band: 1, y: 39, x: 55)>
array([[[     nan,      nan, ...,      nan,      nan],
        [     nan,      nan, ...,      nan,      nan],
        ...,
        [     nan, 0.672842, ...,      nan,      nan],
        [     nan,      nan, ...,      nan,      nan]]])
Coordinates:
  * band     (band) int32 1
  * y        (y) float64 28.75 29.28 29.8 30.33 30.86 ... 47.17 47.7 48.22 48.75
  * x        (x) float64 -100.2 -99.75 -99.25 -98.75 ... -74.25 -73.75 -73.25
Attributes:
    transform:      (0.5, 0.0, -100.5, 0.0, 0.5263157894736842, 28.4868421052...
    crs:            +init=epsg:4326
    res:            (0.5, -0.5263157894736842)
    is_tiled:       0
    nodatavals:     (nan,)
    scales:         (1.0,)
    offsets:        (0.0,)
    descriptions:   (\'yield-soy-noirr\',)
    AREA_OR_POINT:  Area
    grid_mapping:   spatial_ref

标签： python r netcdf python-xarray terra

【解决方案1】：

您可能不再等待任何人在这里回答，但我想在此提供一个答案尝试：

使用terra 读取netCDF 文件时，有几个问题很明显。没有提供坐标参考系统，没有定义范围，并且您的分辨率在 x 和 y 方向上有所不同。就像您在问题中所说的那样，只是想确认一下。

nc <- terra::rast("hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.nc")
nc
#> Error in R_nc4_open: Invalid argument
#> Warning: [rast] GDAL did not find an extent. Cells not equally spaced?
#> class       : SpatRaster 
#> dimensions  : 39, 55, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01818182, 0.02564103  (x, y)
#> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source      : hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.nc:yield-soy-noirr 
#> varname     : yield-soy-noirr 
#> name        : yield-soy-noirr

terra 会引发警告，而不是错误，但仍会创建一个 SpatRaster 对象。

由于我们知道您的光栅尺寸（x：55，y：39）和分辨率（0.5°），我怀疑 - 请证明我错了！ - 这与所提供的范围不一致：

# x
seq(from = -100.25, by = 0.5, length.out = 55) |> range()
#> [1] -100.25  -73.25

# y
seq(from = 28.75, by = 0.5, length.out = 39) |> range()
#> [1] 28.75 47.75

此时你的ymax好像是错误的，但是手动设置extent的时候，得到的分辨率还是0.4909091, 0.4871795 (x, y)，说明extent还是不匹配。玩弄上面的序列，我将length.out 增加到n+1 以便再试一次：

# x
seq(from = -100.25, by = 0.5, length.out = 55+1) |> range()
#> [1] -100.25  -72.75

# y
seq(from = 28.75, by = 0.5, length.out = 39+1) |> range()
#> [1] 28.75 48.25

...手动调整属性后的结果看起来很有希望：

terra::ext(nc) <- c(-100.25, -72.75, 28.75, 48.25)
terra::crs(nc) <- "epsg:4326"

nc
#> class       : SpatRaster 
#> dimensions  : 39, 55, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -100.25, -72.75, 28.75, 48.25  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.nc:yield-soy-noirr 
#> varname     : yield-soy-noirr 
#> name        : yield-soy-noirr

生成的 SpatRaster 对象现在可以写入磁盘：

terra::writeRaster(nc, "hybrid_gfdl-esm4_ssp126_default_yield_soybean_shift_2017-2044.tif")

【讨论】：