在Python中使用rioxarray / rasterio制作空/模板栅格

问题描述 投票:0回答:1

我目前正在尝试将 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)
python r python-xarray terra rasterio
1个回答
0
投票

您的部分困惑在于混淆了 rasterio 和 Rioxarray 的使用。尽管 rioxarray 扩展了 xarray 以使用 rasterio 功能,但最好使用 rasterio 或 rioxarray 打开内容,以执行特定过程。 实际上,我在许多工作流程中都做了类似的事情。

使用rioxarray打开
    “掩模光栅”。您可以轻松地使用其他一些方法来生成“模板”数组。但将其设为 xarray 数据集。 ...
  1. 打开并格式化掩模光栅 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"})

    由于我经常使用数据密集的输入栅格作为我的“模板”数组/数据集,因此我将通过首先剥离数据集变量中的实际值来简化它。
    ...
  2. # 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 数据集:
    打开目标/分析栅格
  1. 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
  2. 执行对齐: ...
  3. grid_ds = grid_ds.rio.reproject_match(mask_ds)

    最后,使用特定类型的 xarray 连接功能,称为 
    combine_by_coords
  4. ,这样我实际上将所有分析数组以及掩码数组作为数据集变量保存在单个 xarray 数据集中。这确保了对齐和随后的轻松处理:
  5. 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)

    	
最新问题
© www.soinside.com 2019 - 2025. All rights reserved.