前言

上一篇文章實作了MODIS Swath資料的重投影，一般來說，在完成重投影之后的任務，要么是對多幅影像的拼接，要么是對多幅影像做均值運算，但這兩個任務，在一定程度上是可以等同的，均值運算，在不考慮質量控制等因素的前提下，其實就是先把所有影像包含的最大地理范圍算出來，然后再在對應的地理位置上，把屬于這一個地理位置上的所有像元給加進來，并統計有效的累加次數，最后求平均即可，對于某些位置上只出現過一次的像元，其實就是把值填進去了，最后累加的次數為1，均值就是像元值本身，而拼接，實際上是相同的道理，比如兩幅有重合區域的影像，重合的部分有多種處理策略，你可以對重合的部分算均值，也可以用下一幅的結果覆寫上一幅的結果，對于不重合的區域，那也不存在處理策略的問題，直接把值填進去即可，這就是為什么我說拼接和均值運算在一定程度上可以等同的原因，
這篇文章里給的代碼，可以根據實際需要使用，拼接和求均值都可以用這個代碼完成，當然有一些使用的前提，后面再詳細說，

Python實作

這個代碼用起來也比較粗暴，整體的邏輯，是從給定的輸入路徑中搜索geotiff檔案（即已經重投影好的結果，是不是MODIS的結果其實無所謂），并給定檔案名中能夠代表日期的字串起止位置，之后程式會對檔案串列中包含的所有日期進行統計，統計完成后按照日期進行回圈處理，最終實作對應日期的資料拼接（均值計算），可以看到這個程式里用了兩個比較關鍵的陣列，一個是data_box_geo_sum，另一個是data_box_geo_num，這兩個陣列是分別用來存盤累加結果和累計有效頻次用的，最終兩個陣列相除就是拼接（均值）結果，以上這些就是這個程式中mosaic_unified_resolution這個函式的功能，具體輸入引數及解釋如下：
input_directory：輸入路徑名稱，路徑下包含了所有待拼接/求均值的geotiff檔案，但記得給定時，不要省略路徑最后的’/’，
output_directory：輸出路徑名稱，不需要自己創建，程式中會檢測路徑是否存在并自動創建，但記得給定時，不要省略路徑最后的’/’，
file_prefix：geotiff檔案的后綴，一般就是’.tiff’或者’.tif’，當然你也可以輸入你自己的后綴，只要你確定它確實是那個后綴并且是geotiff檔案就行，
doy_start_pos：檔案名中能夠代表日期的字串開始位置，
doy_end_pos：檔案名中能夠代表日期的字串結束位置，
注意doy_start_pos和doy_end_pos的輸入是有講究的，對于大部分衛星資料來講，他的檔案名里包含的日期資訊一般都不止一個，拿任意一個MODIS的氣溶膠產品檔案舉例，如MYD04_3K.A2018121.0545.061.2018121172155.hdf，這個檔案里日期資訊就有兩個A2018121.0545是他的觀測起始時間，即2018年的第121天，UTC時間5點45分，2018121172155則是這個產品的生成時間，拆開來就是2018年的第121天，UTC時間17點21分55秒，從這個例子就能看出，如果只是簡單地截取2018121這個字串，給這個函式一個起始位置10結束位置17，那很有可能會把某些有相同特征日期字串的資料拿進來一起拼接了，因為我這個程式只是簡單地在檔案名中尋找是否有完整的目標字串，如果有，那就都作為待處理的目標日期加入file_doy這個串列中，如果doy_start_pos和doy_end_pos選取不當，那可能會加入一些不想要的資料進來，比如這樣一種情況：極有可能某個2018120的觀測結果，是在2018121這天進行的處理，對于這個程式也會把這一天的檔案加進來，如果對天做拼接，那時間上顯然就錯了，所以可以看到我的main()程式里，特意給的起始位置是9，也就是把“A”這個字符加進來了，這樣能避免和后面日期的重復，所以這個doy_start_pos和doy_end_pos怎么取，需要根據實際情況來定，如果現在我們想算2018年所有資料的均值，那可能只需要給個起始位置9結束位置14就行了，這樣程式就會把包含有A2018這個字串的所有資料拿進來進行計算，

import os
import gdal
import math
import numpy


def mosaic_unified_resolution(input_directory, output_directory, file_prefix, doy_start_pos, doy_end_pos):
    if not os.path.exists(output_directory):
        os.mkdir(output_directory)

    file_doy = []
    file_all = os.listdir(input_directory)
    for file_i in file_all:
        doy_temp = file_i[doy_start_pos:doy_end_pos]
        if file_i.endswith(file_prefix) and doy_temp not in file_doy:
            file_doy.append(doy_temp)

    for doy_i in file_doy:
        file_list = []
        for file_i in file_all:
            if file_i.find(doy_i) != -1:
                file_list.append(file_i)
        lon_min = 9999.0
        lon_max = -9999.0
        lat_min = 9999.0
        lat_max = -9999.0
        for file_list_i in file_list:
            data_temp = gdal.Open(input_directory + file_list_i)
            data_col = data_temp.RasterXSize
            data_line = data_temp.RasterYSize
            geo_info = data_temp.GetGeoTransform()
            lon_min_temp = geo_info[0]
            lon_max_temp = lon_min_temp + float(data_col) * geo_info[1]
            lat_max_temp = geo_info[3]
            lat_min_temp = lat_max_temp + float(data_line) * geo_info[5]
            if lon_min_temp < lon_min:
                lon_min = lon_min_temp
            if lon_max_temp > lon_max:
                lon_max = lon_max_temp
            if lat_min_temp < lat_min:
                lat_min = lat_min_temp
            if lat_max_temp > lat_max:
                lat_max = lat_max_temp
        output_resolution_x = geo_info[1]
        output_resolution_y = geo_info[5]
        data_box_geo_col = math.ceil((lon_max - lon_min) / abs(output_resolution_x))
        data_box_geo_line = math.ceil((lat_max - lat_min) / abs(output_resolution_y))
        data_box_geo_sum = numpy.full((data_box_geo_line, data_box_geo_col), 0.0)
        data_box_geo_num = numpy.full((data_box_geo_line, data_box_geo_col), 0.0)
        for file_list_i in file_list:
            data_temp = gdal.Open(input_directory + file_list_i)
            data_geoinfo = data_temp.GetGeoTransform()
            data_array = data_temp.ReadAsArray()
            data_line = data_array.shape[0]
            data_col = data_array.shape[1]
            col_start_pos = int((data_geoinfo[0] - lon_min) / abs(output_resolution_x))
            line_start_pos = int((lat_max - data_geoinfo[3]) / abs(output_resolution_y))
            data_box_geo_sum[line_start_pos:(line_start_pos + data_line),
            col_start_pos:(col_start_pos + data_col)] += data_array
            data_box_geo_num[line_start_pos:(line_start_pos + data_line), col_start_pos:(col_start_pos + data_col)] += (
                    data_array > 0.0)
        data_box_geo_num = (data_box_geo_num > 0.0) * data_box_geo_num + (data_box_geo_num == 0.0) * 1.0
        data_box_geo = data_box_geo_sum / data_box_geo_num
        driver = gdal.GetDriverByName('GTiff')
        out_file = driver.Create(output_directory + doy_i + '.tiff', data_box_geo_col, data_box_geo_line, 1,
                                 gdal.GDT_Float32)
        out_file.GetRasterBand(1).WriteArray(data_box_geo)
        geoinfo_new = [lon_min, output_resolution_x, 0, lat_max, 0, output_resolution_y]
        out_file.SetGeoTransform(geoinfo_new)
        out_file.SetProjection(data_temp.GetProjection())
        print('The mosaic of {} is complete.'.format(doy_i))


def main():
    input_directory = 'O:/coarse_data/chapter_1/MODIS_2018_mod04_3k/geo_out/'
    output_directory = 'O:/coarse_data/chapter_1/MODIS_2018_mod04_3k/geo_out/mosaic/'
    file_prefix = '.tiff'
    doy_start_pos = 9
    doy_end_pos = 17
    mosaic_unified_resolution(input_directory, output_directory, file_prefix, doy_start_pos, doy_end_pos)


if __name__ == '__main__':
    main()

接下來就是程式運行的一些效果，這個是geo_out檔案夾下待拼接的檔案：
這個是運行程式后mosaic檔案夾下生成的結果：
這個是我把截取位置變成9、14之后，算的2018年所有結果的均值：

后記

這個程式要正確運行，numpy和gdal庫必須要有，另外對輸入資料是有要求的：（1）輸入資料必須是已經做過幾何校正或重投影的GeoTIFF資料；（2）每個輸入資料的投影和空間解析度要求統一，比如都是經緯度網格、空間解析度為0.03°；（3）資料中的無效值只有一個，比如0，而不同時有0和-9999，當然具體的無效數值可以自己直接修改，對應像元有效頻次data_box_geo_num計算時data_array > 0.0這個地方，
再強調下這個程式的邏輯是把兩幅影像重合的地方進行均值計算，并作為最終拼接的結果，所以拼接和均值計算都可以用這個程式來實作，如果拼接時需要的是下一幅對上一幅的覆寫，那對應的在data_box_geo_sum這個地方就不能做累加，只能直接等于，可能還需要加一些對日期的判斷，另外data_box_geo_num也就不需要了，整個程式我覺得還是有可讀性，大家可以根據自己的需要進行修改，
另外注意到geo_info[5]這個東西，我這套代碼測驗的時候，用的是上一篇文章里處理出來的幾何校正結果，也就是用GDAL做的，之后再讀出來的地理資訊，geo_info[5]這個實際上代表的是y方向上的解析度，但如果print出來就會發現它是個負值，所以最終歸算最小緯度的時候用的是lat_max_temp + float(data_line) * geo_info[5]，這個和我寫IDL的時候有些不一樣，我不清楚是不是所有GeoTIFF檔案用GDAL讀出來都這樣，
其實整個代碼有很重的個人寫IDL的風格，畢竟才接觸Python不到3天時間，有些習慣很難改過來，