??本文介紹基于Python中的ArcPy模塊，將大量遙感影像文件按照分幅條帶編號(hào)與成像時(shí)間加以分組，并將同一分幅的遙感影像加以每個(gè)8天時(shí)間間隔內(nèi)的鑲嵌拼接的方法。
??首先，來(lái)看一下本文具體的需求。我們現(xiàn)有一個(gè)文件夾，其中含有大量的.tif格式的遙感影像，如下圖所示。首先，每一景遙感影像的文件名中，都有一個(gè)表示成像時(shí)間的字段；例如，下圖中從上往下數(shù)第1景圖像，就是2022年第001天某時(shí)刻的遙感影像，而下圖中從上往下數(shù)第4景圖像，就是2022年第013天某時(shí)刻的遙感影像。
??同時(shí)，這些遙感影像文件的文件名順序還不完全是時(shí)間順序，因?yàn)槠湮募_(kāi)頭還有一些表示其他含義的字段（如傳感器名稱(chēng)），而這些不同字段對(duì)應(yīng)的遙感影像文件同樣具有多個(gè)成像時(shí)間。如下圖所示，可以看到在GF1WFV3傳感器對(duì)應(yīng)的2022年346天遙感影像結(jié)束后，新的GF1WFV4傳感器對(duì)應(yīng)的遙感影像又是從2022年的開(kāi)頭開(kāi)始的?？傊?，就是不能將文件名排序作為遙感影像成像時(shí)間的順序。
??其次，如下圖所示，每一景遙感影像的文件名中還有一個(gè)表示遙感影像分幅的字段；其中，48STA與48STB等都是不同分幅對(duì)應(yīng)的編號(hào)。
??我們希望實(shí)現(xiàn)的是，從2022年第001天開(kāi)始，到第365天結(jié)束，對(duì)于每1個(gè)分幅，將其每1個(gè)8天時(shí)間范圍內(nèi)的所有遙感影像（無(wú)論是來(lái)自哪一個(gè)傳感器）拼接在一起。例如，將分幅為48STA的、成像時(shí)間在001天至008天的遙感影像拼接在一起，然后將009天至016天的拼接在一起，以此類(lèi)推，直到2022年所有分幅為48STA的遙感影像處理完成；隨后再處理48STB的，再以此類(lèi)推，直到全部分幅都處理完成。
??在之前的文章Google Earth Engine谷歌地球引擎求取每隔N天的遙感影像平均數(shù)中，我們介紹過(guò)在GEE中計(jì)算每1個(gè)8天時(shí)間間隔內(nèi)遙感影像數(shù)據(jù)平均值的方法；而這一次我們將基于Python，將每1個(gè)8天時(shí)間間隔內(nèi)遙感影像拼接起來(lái)。
??本文所用到的代碼如下。
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 13 22:44:34 2024
@author: fkxxgis
"""
import re
import os
import arcpy
arcpy.env.workspace = r"F:\Data_Reflectance_Rec\GF\2022"
output_folder = r"F:\Data_Reflectance_Rec\GF\8Days_frame\2022"
image_list = arcpy.ListRasters("*", "tiff")
image_list.sort()
image_dict = {}
for image in image_list:
match = re.search(r"\d{7}", image)
image_date = match.group()
image_year = image_date[0:4]
image_days = image_date[-3:]
match = re.search(r"\.(.{5})\.", image)
image_frame = match.group().strip(".")
dict_idx = (int(image_days) - 1) / 8
dict_key = str(dict_idx) + "/" + image_frame
if dict_key in image_dict:
image_dict[dict_key].append(image)
else:
image_dict[dict_key] = [image]
print image, "is in", dict_key
for dict_idx, image_list_interval in image_dict.items():
dict_idx_split = dict_idx.split("/")
days = int(dict_idx_split[0]) * 8 + 1
frame = dict_idx_split[1]
template_image = image_list_interval[0]
cell_size = arcpy.GetRasterProperties_management(template_image, "CELLSIZEX")
value_type = arcpy.GetRasterProperties_management(template_image, "VALUETYPE")
describe = arcpy.Describe(template_image)
spatial_reference = describe.spatialReference
arcpy.CreateRasterDataset_management(output_folder,
image_year + str(days).zfill(3) + "_" + frame + ".tif",
cell_size.getOutput(0),
"16_BIT_UNSIGNED",
spatial_reference,
4)
print image_year + str(days).zfill(3) + "_" + frame + ".tif", "creation finished."
try:
arcpy.Mosaic_management(image_list_interval,
os.path.join(output_folder, image_year + str(days).zfill(3) + "_" + frame + ".tif"),
"MINIMUM",
background_value = 0,
nodata_value = 0)
# arcpy.MosaicToNewRaster_management(image_list_interval,
# output_folder,
# image_year + str(days).zfill(3) + "_" + frame + ".tif",
# pixel_type = "16_BIT_UNSIGNED",
# number_of_bands = 4,
# mosaic_method = "MINIMUM")
print image_year + str(days).zfill(3) + "_" + frame + ".tif", "mosaic finished."
except arcpy.ExecuteError:
print image_year + str(days).zfill(3) + "_" + frame + ".tif", "had ERROR."
??其中，上述代碼的具體含義如下。
??首先，我們通過(guò)import語(yǔ)句導(dǎo)入所需的模塊。其中，re用于正則表達(dá)式匹配，os用于文件路徑操作，arcpy是ArcGIS的Python模塊，用于處理GIS數(shù)據(jù)。
??隨后，我們通過(guò)env.workspace設(shè)置工作空間，即等待拼接的柵格影像數(shù)據(jù)所在的文件夾路徑；通過(guò)output_folder設(shè)定輸出結(jié)果的文件夾路徑。
??接下來(lái)，基于ListRasters("*", "tiff")獲取待拼接的所有.tif格式柵格文件，并將其排序后存儲(chǔ)在image_list列表中；image_dict是一個(gè)字典，用于存儲(chǔ)柵格影像按日期和分幅號(hào)進(jìn)行分組的結(jié)果，從而將每1種分幅中，處于同1個(gè)8天時(shí)間間隔的遙感影像放在一起；for循環(huán)遍歷image_list中的每個(gè)影像文件，并使用正則表達(dá)式re.search提取影像文件名中的日期信息——其中，需要提取年份image_year和天數(shù)image_days；接下來(lái)，使用正則表達(dá)式re.search提取影像文件名中的分幅號(hào)信息，并根據(jù)天數(shù)和分幅號(hào)生成字典的鍵dict_key；隨后，將影像文件添加到相應(yīng)的字典值中，如果字典鍵已存在，則將影像文件添加到對(duì)應(yīng)的列表中。同時(shí)，打印信息，指示影像文件屬于哪個(gè)字典鍵。
??再次，for循環(huán)遍歷image_dict中的每個(gè)字典鍵和對(duì)應(yīng)的影像文件列表——首先拆分字典鍵，獲取天數(shù)和分幅號(hào)的信息；接下來(lái)，獲取文件列表中第一個(gè)影像文件的信息，如像元大小、值類(lèi)型、空間參考等（因?yàn)楹笃谛枰谄鋪?lái)作為模板圖像）；隨后，使用CreateRasterDataset_management()函數(shù)創(chuàng)建輸出柵格數(shù)據(jù)集，命名規(guī)則為年份+天數(shù)+分幅號(hào)。同時(shí)，打印信息，指示柵格數(shù)據(jù)集創(chuàng)建完成。
??最后，即可使用Mosaic_management()將影像文件列表拼接為一個(gè)柵格數(shù)據(jù)集，命名規(guī)則同上；同時(shí)，打印信息，指示柵格數(shù)據(jù)集拼接完成。如果拼接過(guò)程中出現(xiàn)錯(cuò)誤，則捕獲arcpy.ExecuteError異常，并打印錯(cuò)誤信息。這里之所以需要try和except語(yǔ)句，是因?yàn)橛械?code>8天時(shí)間間隔內(nèi)可能沒(méi)有任何遙感影像數(shù)據(jù)，因此Mosaic_management()函數(shù)可能會(huì)報(bào)錯(cuò)，導(dǎo)致程序終止運(yùn)行。關(guān)于try和except語(yǔ)句的具體用法，我們將在后續(xù)的博客中介紹。
??運(yùn)行上述代碼，首先將看到如下圖所示的界面；表示正在基于遙感影像的文件名，將其放置到不同的字典中——這個(gè)字典就是根據(jù)遙感影像成像時(shí)間與分幅號(hào)來(lái)表示的。
??完成字典的確定后，相同分幅號(hào)且落在同1個(gè)8天時(shí)間間隔內(nèi)的遙感影像數(shù)據(jù)，即可被存入同1個(gè)字典中。接下來(lái)，即可開(kāi)始拼接；如下圖所示。
??完成上述代碼運(yùn)行后，即可在結(jié)果文件夾中看到按照分幅號(hào)與成像時(shí)間拼接好的遙感影像了。因?yàn)槲疫@里當(dāng)初把2022年的拼接結(jié)果誤刪了，所以就截取2021年的數(shù)據(jù)經(jīng)過(guò)上述代碼處理后的結(jié)果，如下圖所示?？梢钥吹剑Y(jié)果已經(jīng)是按照每個(gè)8天的時(shí)間間隔、以及每1種分幅號(hào)拼接好的了。
??至此，大功告成。
