arcpy在arcgis中的求交分析与数据融合

首先这是一个利用arcpy求交、以及数据融合、文件复制操作例子。大概需求是这样的，这里有一个图幅结合表的shapefile面数据，该面属性表数据对应有图幅号，现在利用某一区域（暂且一个县）的面shapefile数据与图幅结合表相交，计算出相交结果的图幅结合表，然后遍历该结果图幅结合表的图幅字段，最后利用利用这些字段信息在指定文件将多个分幅好的DEM数据拼接为该区域（一个县）一大块DEM数据，整个工作流程就结束了。

来看一下相关的文件。如下图所示是一个图幅里面对应的DEM数据文件。

现在利用arcpy的函数Intersect_analysis计算出图幅结合图表相交面，然后利用里面的属性数据找出所有上面所有文件的DEM数据。最后利用MosaicToNewRaster_management函数将上面所有的DEM数据融合到一起。其中利用求交函数得出的DEM文件如下所示。下图的文件夹对应上图中的各个文件。

在MosaicToNewRaster_management函数中传入了，各个文件夹里面的DEM绝对路径的字符串，设置相应的输出路径、文件名、以及合成DEM波段，因为DEM也是可以是一种tif格式数据，这里根据原始数据的设置了32个波段。而如果使用Arcgis的工具箱来合成新的数据，依次选择【Data ManageMent Tools】、【Raster】、【Raster Dataset】、【Mosaic To New Raster】，弹出如下的设置界面对话框。

设置完上面的合成界面后，就可以做数据融合工作。注意上面数据就两个，在这里我们可以将上面选择arcgis的toolbox转为相应的python脚本，供我们使用，可以说非常方便。具体操作如下，操作完上面的融合测试工作后（因为只是选择了两个DEM数据，而我们接下来需要选择一个区域的多个DEM数据），我们在Arcgis依次选择【Geoprocessing】、【Results】

右击刚才的融合操作【Mosaic To New Raster】、选择【Copy As Python Snippet】，就可以将刚才的融合需要输入的参数、输入参数信息以及需要调用的arcpy函数变为python脚本。具体如下图所示。

这也是需要在我们自己脚本中调用的python脚本，经过这样的操作极大方便我们使用arcpy函数。

好了，说了这么多，来看一下我们写的python脚本吧。

# coding:utf-8
import arcpy
import shutil,os


def del_file(path):
    ls = os.listdir(path)
    for i in ls:
        c_path = os.path.join(path, i)
        if os.path.isdir(c_path):
            del_file(c_path)
        else:
            os.remove(c_path)

def list_all_files(rootdir):

       _files = []
       list = os.listdir(rootdir)  # 列出文件夹下所有的目录与文件
       for i in range(0, len(list)):
             path = os.path.join(rootdir, list[i])
             if os.path.isdir(path):
                _files.extend(list_all_files(path))
             if os.path.isfile(path):
                 _files.append(path)
       return _files

def containVarInString(containVar,stringVar):
    try:
        if isinstance(stringVar, str):
            if stringVar.find(containVar):
                return True
            else:
                return False
        else:
            return False
    except Exception,e:
        print e

#http://www.mamicode.com/info-detail-2313739.html
jctb=ur"F:/2018/STTJ/接合图表/gz一万图幅2000.shp"
fwx=ur"F:/2018/STTJ/FWX/qnbuffer/sandubuffer.shp"
resultShp=ur"F:/2018/STTJ/FWX/temp/temp.shp"

resultPath=ur"F:/2018/STTJ/FWX/result/LD/"
dataPath=ur"F:/2018/STTJ/DEM/"

shpRoot=ur"F:/2018/STTJ/FWX/temp/"
#如果之前的结果shp存在,则删除
#if os.path.exists(resultShp):
#    os.remove(resultShp)
#else:

if os.path.exists(shpRoot):
    files= list_all_files(shpRoot)
    for tempfile in  files:
        os.remove(tempfile)


resultShp=ur"F:/2018/STTJ/FWX/temp/temp.shp"


#overlapedShp = arcpy.GetParameterAsText(0)
#resultShp = arcpy.GetParameterAsText(1)

arcpy.Intersect_analysis([jctb,fwx],resultShp,"ALL", "", "INPUT")


theFields = arcpy.ListFields(resultShp)
FieldsArray = []
for Field in theFields:
     FieldsArray.append(Field.aliasName)

#cursor = arcpy.UpdateCursor(resultShp)
cursor = arcpy.da.SearchCursor(resultShp, ['NewMapNo'])

imgArray=[]

for row in cursor:
     dirFileName=row[0]
     tmpFileDir=dataPath+str(dirFileName)

     #如果文件存在
     if os.path.exists(tmpFileDir):
      tmpResultDir=resultPath+dirFileName
      demFileDir =tmpResultDir+"/"+dirFileName+"DEM"
      infoFileDir=tmpResultDir+"/"+"INFO"

      dataResource=tmpFileDir+"/"+dirFileName+"DEM"

      imgArray.append(dataResource);

      #如果结果文件夹不存在，则创建
      if not os.path.exists(tmpResultDir):
          os.makedirs(tmpResultDir)
          os.makedirs(demFileDir)
          os.makedirs(infoFileDir)
          listDirFiles =list_all_files(tmpFileDir)
          flag=0;
          for fileDir in  listDirFiles:

             if ".xls" in str(fileDir):
                 #shutil.copy(fileDir, tmpResultDir)
                 continue
             elif (str(dirFileName)+'DEM'in str(fileDir)) :#DEM目录下的
                 if (".AUX.xml" in str(fileDir)):
                     #shutil.copy(fileDir, tmpResultDir)
                     continue
                 else:
                    #shutil.copy(fileDir,demFileDir)
                    continue
             elif "INFO" in str(fileDir):#INFO目录下的
                 #shutil.copy(fileDir,infoFileDir)
                 continue
             else:#".AUX.xml"
                 #shutil.copy(fileDir, demFileDir)
                 continue

      else:
          #删除之前创建目录下所有的文件
          del_file(tmpResultDir)
          listDirFiles =list_all_files(tmpFileDir)

          t="xx"
          #复制后面的数据
          #shutil.copy(tmpFileDir, tmpResultDir)

strDataRoot=""

for item in  imgArray:
    print item

length=len(imgArray)
for index in range(len(imgArray)):
    if index==0:
     strDataRoot=imgArray[index]
    else:
     strDataRoot = strDataRoot + ";" + imgArray[index]




arcpy.MosaicToNewRaster_management(input_rasters=strDataRoot,
                                   output_location="I:/20180918/New Folder",
                                   raster_dataset_name_with_extension="sandu.tif",
                                   coordinate_system_for_the_raster="#",
                                   pixel_type="32_BIT_FLOAT",cellsize="#",
                                   number_of_bands="1",
                                   mosaic_method="LAST",
                                   mosaic_colormap_mode="FIRST")

至此，整个说明就完了。

---

                                                                                 更多内容，请关注公众号

                                                                          

 

---