Python arcpy创建栅格、批量拼接栅格 - 疯狂学习GIS
source link: https://www.cnblogs.com/fkxxgis/p/17150814.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
Python arcpy创建栅格、批量拼接栅格
本文介绍基于Python语言arcpy
模块,实现栅格影像图层建立与多幅遥感影像数据批量拼接(Mosaic)的操作。
首先,相关操作所需具体代码如下:
import os
import arcpy
file_path="G:/Postgraduate/LAI_Glass_RTlab/A2018161_Dif/DRT/"
out_file_path="G:/Postgraduate/LAI_Glass_RTlab/A2018161_Dif/DRT/"
out_file_name="Global.tif"
file_name_list=os.listdir(file_path)
tif_file_path=file_path+file_name_list[0]
cell_size_x=arcpy.GetRasterProperties_management(tif_file_path,"CELLSIZEX")
cell_size=cell_size_x.getOutput(0)
value_type=arcpy.GetRasterProperties_management(tif_file_path,"VALUETYPE")
describe=arcpy.Describe(tif_file_path)
spatial_reference=describe.spatialReference
arcpy.CreateRasterDataset_management(out_file_path,out_file_name,cell_size,"16_BIT_SIGNED",
spatial_reference,"1")
out_file=out_file_path+out_file_name
for file in file_name_list:
file_path_name=file_path+file
print(file_path_name)
arcpy.Mosaic_management([file_path_name],out_file)
其中,file_path
为存放有多景初始遥感影像的路径格式为.tif
栅格文件(如果不是.tif
格式,例如是.hdf
等文件,需首先进行文件格式的转换);out_file_path
为拼接后所得结果栅格图层的存放路径;out_file_name
为拼接后所得结果栅格图层的文件名称,其可选格式有很多,如下图所示。
在这里,我们默认所得拼接结果图层为一个(也就是file_path
文件夹中全部的待处理遥感影像最终全拼接在一起);如果大家需要使得拼接结果图层是多幅(也就是file_path
文件夹中待处理遥感影像依据区域、时间等分为很多不同的部分,每一部分拼接在一起),可以参考Python GDAL读取栅格数据并基于质量评估波段QA对指定数据加以筛选掩膜,利用其中的循环方式实现需求。
随后,通过os.listdir()
函数获取file_path
路径下的栅格文件,并存储于file_name_list
列表中。
接下来需要创建一个新的栅格图层。之所以要进行这一步骤,是因为本文后期选择用arcpy.Mosaic_management()
函数进行栅格的批量拼接,因此需要首先创建一个新的、空的栅格图层作为拼接的基准。如果大家的需求不是批量拼接栅格数据,而是单纯想利用arcpy
进行新栅格的创建,那就只看这一部分的代码即可。
在这里,我们选择用file_path
路径下的第一个栅格数据(下称“第一栅格”)作为新栅格图层中各项属性(例如像素边长、像素数据格式等)的依据。首先,arcpy.GetRasterProperties_management()
函数获取第一栅格的像素x
边边长;因为一般栅格数据中像素都是正方形,因此我们就通过cell_size=cell_size_x.getOutput(0)
将第一栅格的像素x
边边长作为新栅格图层像素x
边与y
边二者的边长。再利用arcpy.GetRasterProperties_management()
函数获取第一栅格的数据格式;最后利用中间变量describe
获取第一栅格的空间参考信息。
完成以上步骤后,将已获取的第一栅格的各类信息通过函数arcpy.CreateRasterDataset_management()
带入新栅格中。在这里需要注意:尽可能在将要拼接时选择新栅格为"16_BIT_SIGNED"
及以下的数据格式(具体数据格式类别如下图),且将file_path
路径下待拼接的栅格数据的数据格式也全部修改为这一格式;否则可能会由于数据量大而导致拼接过程极慢。我之前就是由于选用了32 bit float
格式的栅格数据进行拼接,导致全球范围的MODIS一个植被产品数据拼接花了将近一天的时间。如果大家的栅格像素数据包含小数,可以通过乘上一个缩放系数的方式进行数据整数化。
代码最后的一个for
循环,就是遍历file_name_list
中的各个栅格数据,并通过arcpy.Mosaic_management()
函数加以拼接即可。
以上,便完成了本次批量拼接的操作。这里还有一点需要注意:由于arcpy
模块的限制,如果大家的Python版本是3.0
及以上,往往不能直接运行上述代码,最好是在ArcMap的Python运行框或其对应IDLE(如下图所示)中运行。
至此,大功告成。
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK