python对ASCLL格式栅格数据进行处理与运算

Question

#利用arcgispro2.5内置的python3.6.9进行如下栅格数据的处理：（或python2.7下的arcpy处理）
根据要求对以下步骤写出详细代码：
1.指定路径1(E:\0104YUZHITEST\0104YUZHITEST)的ASCALL1文件夹下存在以变量I和变量D命名为"I{I}D{D}"的ASCLL格式的栅格文件，其中这两个变量可以分别从路径1下的IALL.xls（列名为qiangdu)和DALL.xls(列名为lishi）表格中读取。每一个栅格文件在arcgispro中打开均是重合的，均有617行，589列，但该ASCLL格式的栅格文件非规则矩形，即存在部分NODATA值，具体的值为-9999. ，有值的每一个栅格都对应一个FS值，类型为小数）
单个栅格文件在arcgispro在打开如下图所示：

在txt记事本中打开如下图所示：

2.由于每一个栅格文件在arcgispro中打开均是重合的，以变量I从小到大的顺序遍历打开指定路径1(E:\0104YUZHITEST)的ASCALL1文件夹下所有不同变量I和变量D组合的栅格文件（一共n个，可以用arcgispro打开），遍历每一个栅格文件下的每一个栅格单元，获得其对应的FS值，为每一个栅格单元创建一个三列的表格（列名分别为M,N,P) ，记录每一个栅格单元记录所打开每一个文件的文件名的变量I值（列名为P）与变量D值（列名为M）以及对应该栅格单元的FS值（列名为N），最终每一个栅格单元的表格都会记录n（n为所打开的栅格文件个数）行数据

3.以上面步骤所获取的每一个栅格单元的记录表中的n个数组（M,N,P)进行以下分析：对于每一个栅格单元，都以P值相等的归为一类，共划分划分为多个(M,N)值坐标序列（即有e个I值，就有e个(M,N)值坐标序列）。分别绘制每个序列的坐标点为折线图（其中M值为横坐标值，N为纵坐标值，并将对应的P值作为对应折线的名称。

4.求解每一个栅格单元上面这些折线与N=1这条直线的交点，横坐标记录为X，与横坐标与之对应的折线名P值记录为纵坐标Y ，由此可获取t个坐标（X,Y），其中t≥0，对每一个栅格单元获取的t个交点（当t≥3时，t＜3则后续直接赋值为0）进行幂函数拟合，获取每一个栅格单元自己的幂函数公式（即I=aD^b)，并计算拟合的相关性系数（R2），所有每一个栅格单元的对应幂函数公式里的参数a、b、拟合的相关性系数R2、参与拟合的点数t，分别以tif格式输出为栅格文件到路径1下
最终达到以下效果：α即a，β即b，该文献链接如下：https://link.springer.com/article/10.1007/s10346-020-01481-9

测试数据可从以下链接获取：链接：https://pan.baidu.com/s/1IuWwAKNg3_SZNVpmTmyDaA?pwd=egop
提取码：egop
--来自百度网盘超级会员V3的分享

梳理一下步骤，请以上面的表述为准：

①首先，您需要指定路径1（E:\0104YUZHITEST）下的ASCALL1文件夹，并确保该文件夹中存在以变量I和变量D命名的ASCLL格式的栅格文件。您可以从路径1下的IALL.xlsx和D.xlsx表格中读取这两个变量的值。

②接下来，您需要使用ArcGIS Pro打开ASCALL1文件夹下的所有包含不同变量I和变量D组合的栅格文件。遍历每个栅格文件的每个栅格单元，获取其对应的FS值，并记录每个栅格单元所打开的文件的变量I值（列名为P）、变量D值（列名为M）以及对应的FS值（列名为N）。每个栅格单元都会记录n个数组（M,N,P），其中n为所打开的栅格文件个数。

③使用上述步骤获取的每个栅格单元的n个数组（M,N,P），进行以下分析：

④对于每个栅格单元，根据不同的P值划分为多个(M,N)值坐标序列。然后，将这些序列的坐标点绘制为折线图，其中M值为横坐标值，N值为纵坐标值，并将对应的P值作为折线的名称。接下来，求解N=1这条直线与各条折线的交点的横坐标X，并记录与横坐标对应的折线名P值作为纵坐标Y。通过这个过程，您可以获取t个坐标点数据（X,Y），其中t≥0。

⑤对于每个栅格单元获取的t个交点，如果t≥3（t＜3则后续直接赋值为0），进行幂函数拟合。通过幂函数拟合，您可以得到每个栅格单元自己的幂函数公式（即I=aD^b），并计算拟合的相关性系数（R2）。

⑥最后，将所有栅格单元的拟合系数a、拟合系数b、拟合相关性系数R2以及参与拟合的点数t保存为四个tif栅格文件，并将其输出到路径1下新建的输出文件目录中。

Answer 1

引自免费微信小程序：皆我百晓生

import arcpy
import os

# 获取变量I和变量D的值

IALL_path = r'E:\0104YUZHITEST\IALL.xlsx'
DALL_path = r'E:\0104YUZHITEST\DALL.xlsx'

I_values = arcpy.da.TableToNumPyArray(IALL_path, "qiangdu")
D_values = arcpy.da.TableToNumPyArray(DALL_path, "lishi")

# 指定路径1下的ASCALL1文件夹，并获取所有栅格文件

input_folder = r'E:\0104YUZHITEST\ASCALL1'
output_folder = r'E:\0104YUZHITEST\output'

# 创建输出目录
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

for i in range(len(I_values)):
    for j in range(len(D_values)):
        I = I_values[i]
        D = D_values[j]
        file_name = f'I{I}D{D}.asc'
        file_path = os.path.join(input_folder, file_name)

        # 打开栅格文件并处理
        try:
            arcpy.env.workspace = output_folder
            arcpy.MakeRasterLayer_management(file_path, "layer")
            arcpy.CopyRaster_management("layer", f"output_{i}_{j}.tif")
            
            # 进行进一步处理和分析
            # ...
            
        except Exception as e:
            print(f"处理文件'{file_name}'时发生错误：{str(e)}")