slope_栅格数据slope_

def main(): print "load modules..." import arcpy import numpy as np # template = r'C:\GeoNet\LoopMultiFolderCalcMeanSameFolder\Test_Data\2004A\{0}162004.tif' template = r'C:\GeoNet\WUE\r{0}_WUE.TIF' nodata = -3.4028235e+38 out_ras = r'C:\GeoNet\WUE\WUE_slope.tif' print "create numpy list..." lst_np_ras = [] for i in range(1, 14): ras_path = template.format("%03d" % (i,)) print " - ", ras_path ras_np = arcpy.RasterToNumPyArray(ras_path) lst_np_ras.append(ras_np) print "read props numpy raster..." ras_np = lst_np_ras[0] rows = ras_np.shape[0] cols = ras_np.shape[1] print "create output numpy array..." ras_path = template.format("%03d" % (1,)) raster = arcpy.Raster(ras_path) ras_np_res = np.ndarray((rows, cols)) # np.ndarray((raster.height, raster.width)) print "loop through pixels..." pix_cnt = 0 for row in range(rows): for col in range(cols): pix_cnt += 1 if pix_cnt % 10000 == 0: print "row:", row, " col:", col, " pixel:", pix_cnt lst_vals = [] for ras_np in lst_np_ras: lst_vals.append(ras_np[row, col]) lst_vals = ReplaceNoData(lst_vals, nodata) # perform calculation on list a, b, sum_abs_dif, avg_abs_dif, cnt_elem = RegressionLine(lst_vals) if cnt_elem > 1: val = a ## print " - ", lst_vals ## print " - slope (a):", a ## print " - intercept (b):", b ## print " - total distance from regression line:", sum_abs_dif ## print " - mean distance from regression line :", avg_abs_dif else: val = nodata try: ras_np_res[row, col] = val except: print "err" pnt = arcpy.Point(raster.extent.XMin, raster.extent.YMin) xcellsize = raster.meanCellWidth ycellsize = raster.meanCellHeight ras_res = arcpy.NumPyArrayToRaster(ras_np_res, lower_left_corner=pnt, x_cell_size=xcellsize, y_cell_size=ycellsize, value_to_nodata=nodata) ras_res.save(out_ras) arcpy.DefineProjection_management(in_dataset=out_ras, coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]") def ReplaceNoData(lst, nodata): res = [] for a in lst: if a < nodata / 2.0: res.append(None) else: res.append(a) return res def RegressionLine(lst): from numpy import arange, array, ones, linalg months = range(1, 14) lst_x = [] lst_y = [] cnt = 0 for a in lst: if not a is None: lst_x.append(months[cnt]) lst_y.append(a) cnt += 1 cnt_elem = len(lst_x) if cnt_elem >= 2: A = array([ lst_x, ones(len(lst_x))]) # linearly generated sequence w = linalg.lstsq(A.T, lst_y)[0] # obtaining the parameters of the regression line a = w[0] # slope b = w[1] # offset # calculate acumulated abs difference from trend line zip_lst = zip(lst_x, lst_y) sum_abs_dif = 0 for xy in zip_lst: x = xy[0] y = xy[1] calc = x * a + b abs_dif = abs(y - calc) sum_abs_dif += abs_dif avg_abs_dif = sum_abs_dif / cnt_elem return a, b, sum_abs_dif, avg_abs_dif, cnt_elem else: return None, None, None, None, None if __name__ == '__main__': main()