如何使用Rasterio更改栅格数据的数据类型

Question

我在使用Python的rasterio包时，遇到了处理没有数据值的问题。具体来说，我想对一个栅格数据应用多边形蒙版。这个特定的栅格是Landsat uint8类型的，有7个波段，并且没有明确定义没有数据值，因为255是没有数据的预留值。但是，有时候uint8数据会从uint16数据压缩而来，255值是有效数据值之一，我不希望将其视为“没有数据”（这些数据属于完整位范围）。如果没有指定该参数，则rasterio的掩码函数默认认为0是“没有数据”值，这种方式同样存在问题，因为0有时被认为是有效数据值。是否有一种方法可以覆盖元数据值的“没有数据”？

我尝试了几种解决方法（详见下文），但都没有成功。

1. 使用rasterio.open（）将uint8数据转换为uint16数据，并将“256”作为没有数据值分配。因为它超出了任何uint8数据的范围，但在uint16数据范围内可接受。这是某些软件程序（如ArcMap）有时处理分配无数据值的方式。 2. 类似于第1步，但尝试在函数中使用rasterio.open（）打开uint8数据并设置“nodata = np.nan”。收到错误消息：“给定的no data值nan超出了其数据类型的有效范围。”尽管在文档中列出nan作为“nodata”参数的有效输入项。 3. 在使用rasterio.mask（）进行掩码处理期间，指定nodata = nan。收到错误消息“无法将fill_value nan转换为dtype”。

import rasterio
import fiona
import numpy as np

fp_src = ''
fp_dst = ''
shape = ''

# get shapes
with fiona.open(shape, 'r') as shapefile:
    geoms = [feature['geometry'] for feature in shapefile]



# Method Number 1
# ~~~~~~~~~~~~~~~~~~~~~~~~~

# open original raster, copy meta & alter dtype
with rasterio.open(fp_src) as src_dataset:
    kwds = src_dataset.profile
    kwds['dtype'] = 'uint16'
    src_meta = src_dataset.meta

# write a new raster with the copied and altered meta
with rasterio.open(fp_dst, 'w', **kwds) as dst_dataset:
    dst_meta = dst_dataset.meta

src_dataset.close()
dst_dataset.close()

img = rasterio.open(fp_dst)

# mask img and set nodata to 256 (out of the uint8 range)
out_image, out_transform = mask(img, geoms, nodata=256)
# out_image output: values outside of the geoms are 256 & values inside are 0.



# Method Number 2
# ~~~~~~~~~~~~~~~~~~~~~~~~~

# open original raster, copy meta & alter dtype
with rasterio.open(fp_src) as src_dataset:
    kwds = src_dataset.profile
    kwds['nodata'] = np.nan
    kwds['dtype'] = 'uint16'
    src_meta = src_dataset.meta

# write a new raster with the copied and altered meta
with rasterio.open(fp_dst, 'w', **kwds) as dst_dataset:
    dst_meta = dst_dataset.meta

src_dataset.close()
dst_dataset.close()

img = rasterio.open(fp_dst)

# mask img and let the mask function default to the image's newly created nodata (np.nan from inside with rastario.open...)
out_image, out_transform = mask(img, geoms)
# out_image output: nodata value, nan, is beyond the valid range of its data type



# Method Number 3
# ~~~~~~~~~~~~~~~~~~~~~~~~~

# mask img and set nodata to nan
out_image, out_transform = mask(fp_src, geoms, nodata=np.nan)
# out_image output: Cannot convert fill_value nan to dtype.

我希望看到给定多边形外的所有像素转换为“无数据”条目，该条目不一定是有效范围的一部分，以便脚本不会将有效值误认为是无数据。

- Mike B

Answer 1

问题似乎是你没有将任何数据写入你的新的和修改过的栅格图像，只有元数据。对于方法1，你可以尝试：

# open original raster, copy meta & alter dtype
with rasterio.open(fp_src) as src_dataset:
    kwds = src_dataset.profile
    kwds['dtype'] = 'uint16'
    src_meta = src_dataset.meta
    data = src_dataset.read()

# write a new raster with the copied and altered meta
with rasterio.open(fp_dst, 'w', **kwds) as dst_dataset:
    dst_meta = dst_dataset.meta
    dst_dataset.write(data)

# not needed since you are using a 'with' statement
# src_dataset.close()
# dst_dataset.close()

with rasterio.open(fp_dst) as img:
    # mask img and set nodata to 256 (out of the uint8 range)
    out_image, out_transform = mask(img, geoms, nodata=256)
    # out_image output: values outside of the geoms should be 256 
    # & values inside should be their actual value.

另一种可能更简单的方法是在rasterio的mask函数中使用filled=False选项。这样，您可以获得一个掩蔽的numpy数组，其中所有在掩蔽几何体之外的值都被掩蔽。

- Sam

Answer 2

问题在于np.nan是一个无法转换为整数的浮点数。以下是我解决这个问题的方法：

with rasterio.open(fp_src) as src_dataset:
    meta = src_dataset.meta
    meta.update(
        {
            "nodata": np.iinfo(src_dataset.dtypes[0]).max
        }
    )
    data = src_dataset.read()

    with rasterio.open(fp_dst, 'w', **meta) as dst_dataset:
        dst_dataset.write(data)

函数np.iinfo(dtype).max可以找到给定整数类型的最大值。我在一个只有1个波段的数据集上使用了这个方法，所以它应该也适用于你的数据。如果不行，请告诉我。

- TaaTaa