关于#arcpy#的问题,如何解决?

将1991-1999.tif、2000-2008.tif,这两个栅格文件,需要通过RasterToNumPyArray (arcpy)、numpy将这两个栅格数据在时间维度进行拼接拼接,使其成为1991-2008.tif从而得到土地利用变化情况。
我的尝试:

arr_1991 = arcpy.RasterToNumPyArray("F:\study\luoyang text\CLCD_v01_1991_albert_henan.tif")
arr_1991 = arcpy.RasterToNumPyArray("F:\study\luoyang text\CLCD_v01_1992_albert_henan.tif")
arr_1991_1992 = np.concatenate((arr_1991, arr_1992),axis=0)
new_raster = arcpy.NumPyArrayToRaster(arr_1991_1992, arcpy.Point(x, y),
                                      cell_w, cell_h)
new_raster.save("F:\study\luoyang text\or_1991_1992.tif")  # 保存为一个新的栅格数据


但是结果不是我想要的,我只是在垂直方向按列进行了拼接,但我我需要在在空间维度进行拼接。

img

我需要的是得到一个图像。我在数组拼接应该出现了问题。

假定灰度像素值0~255表示黑到白共256级灰度
那么
新像素值=(int)((图1像素值+图2像素值)/2)

可以试试看

fcs = arcpy.ListFeatureClasses(feature_type='Point')
 
for fc in fcs:
 
    print(fc)
    # Process: 空间连接
    fieldmappings = arcpy.FieldMappings()
    fieldmappings.addTable(featurename)
    arcpy.SpatialJoin_analysis(featurename, fc, 'test', "JOIN_ONE_TO_ONE", "KEEP_ALL",fieldmappings, "INTERSECT")
 
    #删除上一次的目标要素featurename
    arcpy.DeleteFeatures_management(featurename)
        
    #拷贝这一次的输出要素,重新设定featurename
    arcpy.CopyFeatures_management('test',featurename)
 
    #获取featurename字段名
    # fieldName = []
    # fields = arcpy.ListFields(featurename)
    # for field in fields:
    #     fieldName.append(field.name)
    # print(fieldName)
 
    # 删除这一次的输出
    arcpy.DeleteFeatures_management('test')
 
    # 删除多余的字段
    arcpy.DeleteField_management(in_table=featurename,drop_field=['INTERSECT','TARGET_FID'])[0]
 
    # 增加计数字段,设定为图层名
    arcpy.AddField_management(featurename,fc,'LONG')[0]
 
    # Process: 计算字段 (计算字段) 
    arcpy.CalculateField_management(in_table=featurename, field=fc, expression="!Join_Count!", expression_type="PYTHON3", code_block="", field_type="TEXT")[0]
 
    # 删除join_count字段
    arcpy.DeleteField_management(in_table=featurename,drop_field=['Join_Count'])[0]
 
     #明确当前处理过程
    print('当前处理完{}图层'.format(fc))

想要使用arcpy拼接栅格数据的话,需要首先创建一个新的、空的栅格图层作为拼接的基准。然后使用arcpy.Mosaic_management()函数进行栅格的批量拼接。另外如果数据量大,拼接速度很慢的话,选择新栅格为"16_BIT_SIGNED"及以下的数据格式

空间维度叠加,是要两个tif重叠在一起么?这个是tif数据的坐标相同才可以

以下答案由GPT-4.5Plus大模型(语音)与博主A努力中.......共同编写:

您好,根据您的描述,您遇到的问题是在栅格数据拼接时,方向出现错误,导致拼接后得到的是行方向上的拼接,而不是空间维度上的拼接。
要解决这个问题,您需要:

  1. 使用RasterToNumPyArray()将两个栅格转为数组,这个步骤您已完成。
  2. 对两个数组进行轴置换,使时间轴(数组第0轴)变为空间轴(数组第2轴)。可以使用np.swapaxes()实现。
  3. 在第2轴上对两个数组进行拼接,使用np.concatenate()。
  4. 将拼接后的数组使用NumPyArrayToRaster()转换回栅格。

具体代码如下:

python
arr_1991 = arcpy.RasterToNumPyArray("1991-1999.tif")
arr_2000 = arcpy.RasterToNumPyArray("2000-2008.tif")

# 轴置换,交换轴0和轴2
arr_1991 = np.swapaxes(arr_1991, 0, 2)  
arr_2000 = np.swapaxes(arr_2000, 0, 2)

# 在轴2上拼接
arr = np.concatenate((arr_1991, arr_2000), axis=2)   

# 转换回栅格       
new_raster = arcpy.NumPyArrayToRaster(arr, arcpy.Point(x, y), 
                                      cell_w, cell_h)
new_raster.save("1991-2008.tif")

这段代码会先将两个栅格数组的时间轴和空间轴交换,然后在空间轴上进行拼接,最终转换成一张包含时间段1991-2008的新栅格。