我正在尝试将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")