1

我制作了一个脚本，该脚本导入了一个栅格（tif），将其转换为一个数组，经过一些计算，将结果保存在一个新的栅格文件中。但是，我面临的问题是，在第 94 列之后，输出将始终相同，因为邻居 [] 保持为空。我检查了 arcmap，第 94 列之后的单元格确实有值。我发现 95 也是行数。似乎我的脚本在创建正方形后停止正常工作。有谁知道如何解决这个问题？我对python还是很陌生！

这是脚本

import arcpy
import numpy
import itertools
from operator import itemgetter
import copy

Input_raster = 'D:\Master Thesis VSpijker\Python_Scripts\Rasters_Python\Clip_100x100.tif'
Output_raster = 'D:\Master Thesis VSpijker\Python_Scripts\Rasters_Python\is100x100_Raster5.tif'

##Import raster using numpy
PyRaster = arcpy.RasterToNumPyArray(Input_raster)
##A copy of the PyRaster is created to replace the geomorphon values in
LandformRaster = copy.deepcopy(PyRaster)


def Get_neighbors(Raster, dist=1):

     for j in range(len(Raster)):
        for i in range(len(Raster[j])):

##          Save coordinates for every run to replace them in new raster
            x = j
            y = i
##          Extract the neighbor values of a raster cell
            neighbors = []
            i_min = max(0, i-dist)
            i_max = i+dist+1
            j_min = max(0, j-dist)
            j_max = j+dist+1
            print i_min, i_max, j_min, j_max
            for row in Raster[i_min:i_max]:
                neighbors.append(row[j_min:j_max])
##          From list of list to single list with target cell value included
            print neighbors    
            OL_NB = list(itertools.chain(*neighbors))

##          Split list in list with neighbors and target cell value
            if len(OL_NB) == 4:
                NB_Values = itemgetter(1,2,3)(OL_NB)
                Cell_Value = itemgetter(0)(OL_NB)
            elif len(OL_NB) == 6 and len(neighbors) == 2:
                NB_Values = itemgetter(0,2,3,4,5)(OL_NB)
                Cell_Value = itemgetter(1)(OL_NB)
            elif len(OL_NB) == 6 and len(neighbors) == 3:
                NB_Values = itemgetter(0,1,3,4,5)(OL_NB)
                Cell_Value = itemgetter(2)(OL_NB)
            elif len(OL_NB) == 9:
                NB_Values = itemgetter(0,1,2,3,5,6,7,8)(OL_NB)
                Cell_Value = itemgetter(4)(OL_NB)
            else:
                0

##          Determine from every neighbor if the value is higher, lower or equal to the target cell value
            pos = 0
            neg = 0
            equal = 0
            for i in NB_Values:
                if i > Cell_Value:
                    pos = pos + 1
                elif i < Cell_Value:
                    neg = neg + 1
                else:
                    equal = equal + 1

##          turn pos and neg values into a tring value wich can be looked up in a dictionary
            LF_Value = [str(neg), str(pos)]

            LF_Value_Str = "".join(map(str, LF_Value))

##          Dictionary with possible landforms and their codes. 
            Landform_Classes = {'00': 2, '01': 2, '02': 2, '10': 2, '11': 2, '20': 2, '03': 12, '04': 12, '12': 12, '13': 12, '14': 12,
                            '05': 5, '06': 5, '07': 5, '15': 5, '16': 5, '17': 5, '26': 5, '08': 7, '21': 11, '30': 11, '31': 11,
                            '40': 11, '41': 11, '22': 1, '23': 1, '32': 1, '33': 1, '34': 1, '43': 1, '44': 1, '24': 9, '25': 9,
                            '35': 9, '42': 10, '52': 10, '53': 10, '50': 6, '51': 6, '60': 6, '61': 6, '62': 6, '70': 6, '71': 6, '80': 8}
##          table wich shows which numeber is which landform

##            |1|     Slope 5-15 degrees
##            |2|     Slope < 5 degrees
##            |3|     Slope 15-45 degrees
##            |4|     Slope > 45 degrees
##            |5|     Valley
##            |6|     Ridge
##            |7|     Pit
##            |8|     Peak
##            |9|     Hollow
##            |10|    Spur
##            |11|    Shoulder
##            |12|    Footslope


            landform_value = Landform_Classes[LF_Value_Str]


            LandformRaster[x, y] = landform_value


##  Define the projection system and save the new raster in order to import to arcmap again

dsc=arcpy.Describe(Input_raster)
sr=dsc.SpatialReference
ext=dsc.Extent
ll=arcpy.Point(ext.XMin,ext.YMin)

is100x100_Raster = arcpy.NumPyArrayToRaster(LandformRaster,ll,dsc.meanCellWidth,dsc.meanCellHeight)
arcpy.DefineProjection_management(is100x100_Raster, sr)

is100x100_Raster.save(Output_raster)

已经非常感谢了

文森特

4

0 回答 0