我目前正在尝试将 R 代码转换为 Python。我正在尝试制作一个“模板”栅格文件,然后我将在我的工作流程中“重新投影_匹配”许多其他栅格,以便它们全部对齐。
我可以使用 terra 包在 R 中执行此操作,但不确定如何使用 roxarray 包在 python 中执行此操作。我希望有人能够提供帮助?这是我的 R 代码:
# load lib
library("terra")
# Make an empty raster. By default is has extent of the earth and projection wgs84. Setting the resolution means the rows and cols are automatically calculated.
a <- rast(res = 1)
# Fill it with some random integer values
vals(a) <- as.integer(runif(nrow(a) * ncol(a), 1, 10))
# Write it out to file as an
writeRaster(a, "my_integer_rast.tif", wopt = list(datatype = "INT1U"))
我已经走到了这一步,然后有点迷失/困惑,并认为我应该寻求帮助:
import xarray
import rioxarray
data = [[1, 2, 3], [4, 5, 6], [7, 8, 9]] # Example data
coords = {'x': [0, target_res, target_res*2], 'y': [0, target_res, target_res*2]}
ds = xr.Dataset({'data': (['y', 'x'], data)}, coords=coords)
ds.rio.set_spatial_dims(x_dim='x', y_dim='y', inplace=True)
ds.rio.write_crs("EPSG:4326", inplace=True)
ds.rio.write_transform(inplace=True)
with rasterio.open(
"test.tif, 'w',
driver='GTiff',
transform = ds.rio.transform(),
crs=ds.rio.crs,
dtype=rasterio.float32,
#nodata=ds.rio.nodata,
count=1,
width=ref.rio.width,
height=ref.rio.height) as dst:
dst.write(ds.values)
您的部分困惑在于混淆了 rasterio 和 Rioxarray 的使用。尽管 rioxarray 扩展了 xarray 以使用 rasterio 功能,但最好使用 rasterio 或 rioxarray 打开内容,以执行特定过程。 实际上,我在许多工作流程中都做了类似的事情。
使用rioxarray打开
打开并格式化掩模光栅
print('Starting mask')
mask_array = rioxr.open_rasterio(self.masking_raster, band_as_variable=True, masked=True, chunks=self.chunk_return)
if not mask_array.rio.crs:
mask_array.rio.set_crs(self.crs, inplace=True)
mask_array.rio.write_crs(self.crs, inplace=True)
mask_crs = mask_array.rio.crs
print(f"Mask CRS: {mask_crs}")
mask_bounds = mask_array.rio.bounds()
mask_ds = mask_array.rename_vars({"band_1": "mask_values"})
由于我经常使用数据密集的输入栅格作为我的“模板”数组/数据集,因此我将通过首先剥离数据集变量中的实际值来简化它。 ...
# Simplify mask dataset
if mask_ds.nbytes / 1e+9 > 90:
mask_ds = self.mask_ds_values(mask_ds, "mask_values")
print(f'Mask: \n----\n{mask_ds}')
然后,打开栅格以“执行工作”,也作为 xarray 数据集:
grid_array = rioxr.open_rasterio(self.raw_raster, band_as_variable=True)
grid_array.rio.write_crs(self.crs, inplace=True)
grid_ds = grid_array.rename_vars({"band_1": "wse"})
if not grid_ds.rio.crs:
grid_ds.rio.set_crs(self.crs, inplace=True)
使用reproject.mask
grid_ds = grid_ds.rio.reproject_match(mask_ds)
最后,使用特定类型的 xarray 连接功能,称为combine_by_coords
to_combine = to_combine + [grid_ds, mask_ds]
comb_ds = xr.combine_by_coords(to_combine, compat="no_conflicts",
combine_attrs="drop_conflicts")
comb_ds = comb_ds.chunk(self.chunk_return)