有没有办法将 numpy.where() 用于将 NaN 值作为无数据的栅格数据？答案

【问题标题】：is there a way to use numpy.where() for a raster data with NaN values as no-data?有没有办法将 numpy.where() 用于将 NaN 值作为无数据的栅格数据？
【发布时间】：2020-08-09 20:27:33
【问题描述】：

我有一个包含 NaN 值作为无数据的栅格数据。我想从中计算新的栅格，例如 if raster==0 do statement1，if raster==1 do statement2，如果 raster 在 0 和 1 之间做 statement3，否则不要更改值。如何使用 numpy.where() 函数来做到这一点？

这是我的代码：

import os
import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

os.listdir('../NDVI_output')

ndvi1 = rasterio.open("../NDVI_output/NDVI.tiff")

min_value = ndvi_s = np.nanmin(ndvi) #NDVI of Bare soil
max_value = ndvi_v = np.nanmax(ndvi) #NDVI of full vegetation cover

fvc = (ndvi-ndvi_s)/(ndvi_v-ndvi_s) #fvc: Fractional Vegetation Cover

band4 = rasterio.open('../TOAreflectance_output/TOAref_B4.tiff')
toaRef_red = band4.read(1).astype('float64')
emiss = np.where((fvc == 1.).any(), 0.99,
                 (np.where((fvc == 0.).any(), 0.979-0.046*toaRef_red,
                           (np.where((0.<fvc<1.).any(), 0.971*(1-fvc)+0.987*fvc, fvc)))))

【问题讨论】：

标签： python python-3.x numpy nan rasterio

【解决方案1】：

如果raster 是一个数组，

raster == x 给出一个布尔掩码，形状与raster 相同，指示raster 的哪些元素（在您的情况下为像素）等于x
np.where(arr) 给出数组 arr 的元素索引，其计算结果为真。因此，np.where(raster == x) 给出了raster 中等于x 的像素索引。
当且仅当arr 的至少一个元素评估为真时，np.any(arr) 返回 True。因此，np.any(raster == x) 会告诉您raster 的至少一个像素是否为 x。

假设fvc 和toaRef_red 具有相同的形状，并且您想创建一个新数组emiss 用于发射，如果fvc 为1，则将其设置为0.99，如果fvc 为0，则设置为0.979 - 0.046 * toaRef_red , 到 0.971 * (1 - fvc) + 0.987 * fvc 如果 0 fvc

emiss = np.full(ndvi.shape, np.nan)    # create new array filled with nan
emiss[fvc == 1] = 0.99
mask = fvc == 0
emiss[mask] = 0.979 - 0.046 * toaRef_red[mask]
mask = (fvc > 0) & (fvc < 1)
emiss[mask] = 0.971 * (1 - fvc[mask]) + 0.987 * fvc[mask]

这与：

emiss = np.full(ndvi.shape, np.nan)    # create new array filled with nan
emiss[np.where(fvc == 1)] = 0.99
idx = np.where(fvc == 0)
emiss[idx] = 0.979 - 0.046 * toaRef_red[idx]
idx = np.where((fvc > 0) & (fvc < 1))
emiss[idx] = 0.971 * (1 - fvc[idx]) + 0.987 * fvc[idx]

显然，后者是多余的。这里不需要np.where。

【讨论】：