ラスター (tif) をインポートして配列に変換し、いくつかの計算を行った後、結果を新しいラスター ファイルに保存するスクリプトを作成しました。ただし、私が直面している問題は、列 94 の後、neighbor[] が空のままであるため、出力が常に同じになることです。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)
もうどうもありがとう
ヴィンセント