Python遥感开发之批量镶嵌
- 1.ArcGis镶嵌
- 2.Arcpy实现镶嵌
- 2.1 Arcpy实现单个镶嵌
- 2.2 Arcpy实现批量镶嵌
- 3.GDAL实现镶嵌
前言:主要介绍了遥感数据的镶嵌,其中包括使用ArcGis如何完成镶嵌,如何使用Arcpy和GDAL完成镶嵌。
1.ArcGis镶嵌
是ArcGis中的非常重要的一个功能,具体Arcgis中的使用可以参考以下博客。
《Arcgis使用教程(五)ARCGIS空间数据处理之影像镶嵌(拼接)与裁剪_arcgis镶嵌工具在哪里-CSDN博客》
简单理解就是由多个TIF合成一个TIF,再合成一个TIF的过程中,可以选择最大值法、中值法、均值法、最小值法等。至于为什么要进行合成,多数是因为单幅影像缺失比较大,需要多幅影像合成才可以,如下图所示!
2.Arcpy实现镶嵌
2.1 Arcpy实现单个镶嵌
其中32_BIT_FLOAT表示浮点型,16_BIT_UNSIGNED表示整型,MAXIMUM表示最大值合成,MEAN表示均值合成。
# encoding:utf-8
import os
from arcpy.management import MosaicToNewRaster
if __name__ == "__main__":
file = unicode(r"E:\AAWORK\work\b1\20220816", "utf-8")
out = unicode(r"E:\AAWORK\work\202208合成\bb16arcpy", "utf-8")
os.chdir(file)
names = os.listdir(file)
list1 = []
for name in names:
list1.append(file+"\\"+name)
print list1
# mosaic = MosaicToNewRaster(list1, out, "202208_b20.tif", pixel_type="32_BIT_FLOAT",number_of_bands=1, mosaic_method="MAXIMUM")
# mosaic = MosaicToNewRaster(list1, out, "20220816_b4.tif", pixel_type="16_BIT_UNSIGNED",number_of_bands=1, mosaic_method="MAXIMUM")
mosaic = MosaicToNewRaster(list1, out, "20220816_mean_b1.tif", pixel_type="16_BIT_UNSIGNED",number_of_bands=1, mosaic_method="MEAN")
2.2 Arcpy实现批量镶嵌
在2.1的基础上增加了一个for遍历实现批量处理
# encoding:utf-8
import os
from arcpy.management import MosaicToNewRaster
if __name__ == "__main__":
name_time = '202108'
list_band_name = ["b1","b2","b3","b4","b6","b7","b19"]
out = r"E:\AAWORK\work\{0}合成\bb31mean".format(name_time)
out = unicode(out, "utf-8")
for band_name in list_band_name:
file = r"E:\AAWORK\work\{0}\{1}".format(band_name,name_time)
file = unicode(file, "utf-8")
os.chdir(file)
names = os.listdir(file)
list1 = []
for name in names:
list1.append(file+"\\"+name)
file_name = "{0}_mean_{1}.tif".format(name_time,band_name)
mosaic = MosaicToNewRaster(list1, out, file_name, pixel_type="16_BIT_UNSIGNED",number_of_bands=1, mosaic_method="MEAN")
# mosaic = MosaicToNewRaster(list1, out, file_name, pixel_type="32_BIT_FLOAT",number_of_bands=1, mosaic_method="MEAN")
print file_name
3.GDAL实现镶嵌
借助GDAL中的ReprojectImage方法实现镶嵌
import os
from osgeo import gdal
from osgeo import gdalconst
def get_data_list(file_path, out=""):
list1 = [] # 文件的完整路径
if os.path.isdir(file_path):
fileList = os.listdir(file_path)
if out != "":
for f in fileList:
out_data = out + "\\" + f
out_data = out_data.replace(".HDF", "_ndvi.tif")
list1.append(out_data)
else:
for f in fileList:
pre_data = file_path + '\\' + f # 文件的完整路径
list1.append(pre_data)
return list1
if __name__ == '__main__':
inputfile_path = r"E:\AAWORK\work\b19\202208"
output_tif = r"E:\AAWORK\work\202208合成\bb16\20220801-16_b19.tif"
input_tiffs = get_data_list(inputfile_path)
# 打开第一个TIFF文件以获取参考信息
ref_ds = gdal.Open(input_tiffs[0], gdalconst.GA_ReadOnly)
driver = gdal.GetDriverByName("GTiff")
# 创建输出TIFF文件
output_ds = driver.Create(output_tif, ref_ds.RasterXSize, ref_ds.RasterYSize, 1, gdalconst.GDT_Int32)
output_ds.SetProjection(ref_ds.GetProjection())
output_ds.SetGeoTransform(ref_ds.GetGeoTransform())
# 逐一读取和镶嵌TIFF文件
for input_tiff in input_tiffs:
input_ds = gdal.Open(input_tiff, gdalconst.GA_ReadOnly)
band = input_ds.GetRasterBand(1)
# 镶嵌到输出文件并采用最大值法
gdal.ReprojectImage(input_ds, output_ds, None, None, gdalconst.GRA_Max)
# 释放内存
input_ds = None
# 关闭输出文件
output_ds = None
print("TIFF文件已镶嵌到 %s,采用最大值法。" % output_tif)