我对python和numpy非常陌生,如果我误用了一些术语,我很抱歉。
我已经将光栅转换成二维数字阵列,希望能快速高效地对其进行计算。
试图通过实例澄清:
array([[ 4, 1 , 3 , 2] dtype=float32)我希望输出值(只需手工进行第一行)读到:
array([[ 10, 1 , 6 , 3], etc.对怎么做有什么想法吗?
提前感谢!
对于任何感兴趣的人来说,即将完成的脚本:
#Generate Cumulative Thresholds
#5/15/14
import os
import sys
import arcpy
import numpy as np
#Enable overwriting output data
arcpy.env.overwriteOutput=True
#Set working directory
os.chdir("E:/NSF Project/Salamander_Data/Continuous_Rasters/Canadian_GCM/2020/A2A/")
#Set geoprocessing variables
inRaster = "zero_eurycea_cirrigera_CA2A2020.tif"
des = arcpy.Describe(inRaster)
sr = des.SpatialReference
ext = des.Extent
ll = arcpy.Point(ext.XMin,ext.YMin)
#Convert GeoTIFF to numpy array
a = arcpy.RasterToNumPyArray(inRaster)
#Flatten for calculations
a.flatten()
#Find unique values, and record their indices to a separate object
a_unq, a_inv = np.unique(a, return_inverse=True)
#Count occurences of array indices
a_cnt = np.bincount(a_inv)
#Cumulatively sum the unique values multiplied by the number of
#occurences, arrange sums as initial array
b = np.cumsum(a_unq * a_cnt)[a_inv]
#Divide all values by 10 (reverses earlier multiplication done to
#facilitate accurate translation of ASCII scientific notation
#values < 1 to array)
b /= 10
#Rescale values between 1 and 100
maxval = np.amax(b)
b /= maxval
b *= 100
#Restore flattened array to shape of initial array
c = b.reshape(a.shape)
#Convert the array back to raster format
outRaster = arcpy.NumPyArrayToRaster(c,ll,des.meanCellWidth,des.meanCellHeight)
#Set output projection to match input
arcpy.DefineProjection_management(outRaster, sr)
#Save the raster as a TIFF
outRaster.save("C:/Users/mkcarte2/Desktop/TestData/outRaster.tif")
sys.exit()发布于 2014-05-14 20:02:31
根据您想要处理重复的方式,这可能会起作用:
In [40]: a
Out[40]: array([4, 4, 2, 1, 0, 3, 3, 1, 0, 2])
In [41]: a_unq, a_inv = np.unique(a, return_inverse=True)
In [42]: a_cnt = np.bincount(a_inv)
In [44]: np.cumsum(a_unq * a_cnt)[a_inv]
Out[44]: array([20, 20, 6, 2, 0, 12, 12, 2, 0, 6], dtype=int64)当然,a是您的数组扁平化的地方,然后您必须将其重塑为原始形状。
当然,一旦numpy 1.9退出,您就可以将上面的41和42线压缩成一个更快的单行:
a_unq, a_inv, a_cnt = np.unique(a, return_inverse=True, return_counts=True)发布于 2014-05-14 19:10:09
编辑:
这很难看,但我认为它终于奏效了:
import numpy as np
def cond_cum_sum(my_array):
my_list = []
prev = -np.inf
prev_sum = 0
for ele in my_array:
if prev != ele:
prev_sum += ele
my_list.append(prev_sum)
prev = ele
return np.array(my_list)
a = np.array([[4,2,2,3],
[9,0,5,2]], dtype=np.float32)
flat_a = a.flatten()
flat_a.sort()
temp = np.argsort(a.ravel())
cum_sums = cond_cum_sum(flat_a)
result_1 = np.zeros(len(flat_a))
result_1[temp] = cum_sums
result = result_1.reshape(a.shape)结果:
>>> result
array([[ 9., 2., 2., 5.],
[ 23., 0., 14., 2.]])发布于 2014-05-14 20:16:12
使用更少的矮小和更多的python:
a = np.array([[4,2,2,3],
[9,0,5,2]], dtype=np.float32)
np.array([[sum(x for x in arr if x <= subarr) for subarr in arr] for arr in a])
# array([[ 11., 4., 4., 7.],
# [ 16., 0., 7., 2.]])如果和只考虑项目一次,无论项目出现多少,那么,
np.array([[sum(set(x for x in arr if x <= subarr)) for subarr in arr] for arr in a])
# array([[ 9., 2., 2., 5.],
# [ 16., 0., 7., 2.]])https://stackoverflow.com/questions/23662704
复制相似问题