从netcdf到光栅(tif)的转换会改变分辨率和范围吗



我正在尝试将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中使用不同的光栅包,但它们返回时会显示一条消息,说单元格的间距不相等。这很可能是文件(附在这里(的一个问题,但我不知道在哪里以及如何处理。

复制代码:


# 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

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

当使用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 + 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")

最新更新