有没有办法将 numpy.where() 用于 NaN 值作为无数据的栅格数据?
is there a way to use numpy.where() for a raster data with NaN values as no-data?
我有一个栅格数据,其中包含作为无数据的 NaN 值。我想从中计算新栅格,如果栅格==0 做语句 1,如果栅格==1 做语句 2,如果栅格在 0 和 1 之间做语句 3,否则不要更改值。如何使用 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)))))
如果raster
是一个数组,
raster == x
给出一个布尔掩码,其形状与 raster
相同,指示 raster
的哪些元素(在您的情况下为像素)等于 x
np.where(arr)
给出数组 arr
中计算结果为真的元素的索引。因此,np.where(raster == x)
给出了 raster
中等于 x
. 的像素索引
np.any(arr)
returns 当且仅当 arr
的至少一个元素计算为真时为真。因此,np.any(raster == x)
告诉您 raster
的至少一个像素是否为 x。
假设 fvc
和 toaRef_red
具有相同的形状,并且您想创建一个新数组 emiss
用于发射,如果 fvc
为 1,则将其设置为 0.99,如果 fvc
为 0,则为 0.979 - 0.046 * toaRef_red
,如果 0 < fvc
< 1,则为 0.971 * (1 - fvc) + 0.987 * fvc
,否则为 NaN,您可以执行以下操作:
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
。
我有一个栅格数据,其中包含作为无数据的 NaN 值。我想从中计算新栅格,如果栅格==0 做语句 1,如果栅格==1 做语句 2,如果栅格在 0 和 1 之间做语句 3,否则不要更改值。如何使用 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)))))
如果raster
是一个数组,
raster == x
给出一个布尔掩码,其形状与raster
相同,指示raster
的哪些元素(在您的情况下为像素)等于x
np.where(arr)
给出数组arr
中计算结果为真的元素的索引。因此,np.where(raster == x)
给出了raster
中等于x
. 的像素索引
np.any(arr)
returns 当且仅当arr
的至少一个元素计算为真时为真。因此,np.any(raster == x)
告诉您raster
的至少一个像素是否为 x。
假设 fvc
和 toaRef_red
具有相同的形状,并且您想创建一个新数组 emiss
用于发射,如果 fvc
为 1,则将其设置为 0.99,如果 fvc
为 0,则为 0.979 - 0.046 * toaRef_red
,如果 0 < fvc
< 1,则为 0.971 * (1 - fvc) + 0.987 * fvc
,否则为 NaN,您可以执行以下操作:
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
。