我一年numpy
中的每一天都有 365 个二维数组,显示如下图像:
我将它们全部堆叠在一个 3d numpy 数组中。像素值代表我要去掉的云,我想搜索前 7 天或后 7 天(前 7 层,后 7 层)以找到除云之外的值,然后替换云该像素的其他可能值的值(相应像素在其他天/层中经历的值)。
我是python新手,有点迷茫。
有任何想法吗?
谢谢
我一年numpy
中的每一天都有 365 个二维数组,显示如下图像:
我将它们全部堆叠在一个 3d numpy 数组中。像素值代表我要去掉的云,我想搜索前 7 天或后 7 天(前 7 层,后 7 层)以找到除云之外的值,然后替换云该像素的其他可能值的值(相应像素在其他天/层中经历的值)。
我是python新手,有点迷茫。
有任何想法吗?
谢谢
您实际上是在尝试为您的数组编写一个过滤器。
首先,您需要编写一个函数,当给定一个值数组时,中间一个是当前检查的元素,将返回这些值的一些计算。在您的情况下,该函数将期望采用一维数组并返回最接近非云的中间索引的元素:
import numpy as np
from scipy.ndimage.filters import generic_filter
_cloud = -1
def findNearestNonCloud(elements):
middleIndex = len(elements) / 2
if elements[middleIndex] != _cloud:
return elements[middleIndex] # middle value is not cloud
nonCloudIndices, = np.where(elements != _cloud)
if len(nonCloudIndices) == 0:
return elements[middleIndex] # all values were cloud
prevNonCloudIndex = np.where(nonCloudIndices < middleIndex,
nonCloudIndices, -1).max()
nextNonCloudIndex = -np.where(nonCloudIndices > middleIndex,
-nonCloudIndices, 1).min()
# -1 means no non-cloud index
# pick index closest to middle index
if (abs(prevNonCloudIndex - middleIndex)
<= abs(nextNonCloudIndex - middleIndex)):
return elements[prevNonCloudIndex]
else:
return elements[nextNonCloudIndex]
现在您需要将此函数应用于您感兴趣的元素。为此,您需要一个掩码来指示您对特定元素感兴趣的其他元素。
from scipy.ndimage.filters import generic_filter
# creates 5 days worth of a 3x3 plot of land
input = np.ones((5, 3, 3)) * _cloud
input[0,:,:] = 10 # set first "image" to all be 10s
input[4,0,0] = 12 # uppper left corner of fourth image is set to 12
print "input data\n", input, "\n"
mask = (5, 1, 1)
# mask represents looking at the present day, 2 days in the future and 2 days in
# the past for 5 days in total.
print "result\n", generic_filter(input, findNearestNonCloud, size=mask)
# second and third images should mirror first image,
# except upper left corner of third image should be 12
我通过这个解决了它:
interpdata = []
j = 0
for i in stack:
try:
temp = np.where( stack[j] == 50, stack[j-1], modis[j] )
temp = np.where( temp == 50, stack[j+1], temp )
temp = np.where( temp == 50, stack[j-2], temp )
temp = np.where( temp == 50, stack[j+2], temp )
temp = np.where( temp == 50, stack[j-3], temp )
temp = np.where( temp == 50, stack[j+3], temp )
temp = np.where( temp == 50, stack[j-4], temp )
temp = np.where( temp == 50, stack[j+4], temp )
except IndexError:
print 'IndexError Passed'
pass
else:
pass
interpdata [j, :, :] = temp
j = j + 1
我认为您可以执行以下操作:
data = somehow_get_your_3d_data() #indexed as [day_of_year,y,x]
for i,dat in enumerate(data):
weeks2 = data[max(i-7,i):min(i+7,len(data)), ... ]
new_value = get_new_value(weeks2) #get value from weeks2 here somehow
dat[dat == cloud_value] = new_value