版權聲明：本文為部落客原創文章，遵循 CC 4.0 BY-SA 版權協定，轉載請附上原文出處連結和本聲明。
最近在處理MCD19A2資料，從我上一篇部落格得知，處理之後的MCD19A2資料的存儲方式是4波段的栅格檔案，一個波段對應一個時間，但是檢視每一個時間（波段）的有效值發現并不是能全部分布，如下圖：
第一波段： 第二波段：
第三波段： 第四波段：
一般來說，第二波段的有效值最多，如果嫌麻煩可以直接提取第二波段的資料作為結果。我這裡想把4個波段做一個有效值的均值計算。
比如說這個栅格資料是4個波段，但同一處像元上不是每一個值都是有值的，見下圖：
可以看見第一和第四波段上的值是nodata，是無效值，這兩個值是不參與計算的，也就是說如果要算這個像元的有效值均值，隻參與計算的是第三波段和第二波段的值，即（60+119）/2=89.5，而不是（60+119）/4=44.75
1.python中實作兩景影像的有效值均值計算。
這裡的測試資料我是拿MRT校正出來的MCD19A2資料，B2_TIF和B3_TIF就是校正出來的第二波段和第三波段的TIF，這裡的無效值為-28672。
import numpy as np
from osgeo import gdal
B2_TIF=r'G:\0test\mcd19a2.Optical_Depth_055.Orbits_02.tif'
B3_TIF=r'G:\0test\mcd19a2.Optical_Depth_055.Orbits_03.tif'
B2_data=gdal.Open(B2_TIF)
B3_data=gdal.Open(B3_TIF)
#讀取B2TIF的值
rows1 = B2_data.RasterYSize
cols1 = B2_data.RasterXSize
band1 = B2_data.GetRasterBand(1) # 1100 2305
b2_ds = band1.ReadAsArray(0, 0, cols1, rows1)
#讀取B3TIF的值
rows2 = B3_data.RasterYSize
cols2 = B3_data.RasterXSize
band2 = B3_data.GetRasterBand(1) # 1100 2305
b3_ds = band2.ReadAsArray(0, 0, cols1, rows1)
out_ds=np.zeros((rows1,cols1))
#因為我這裡的無效值全都是小于0的值，是以我根據0值這個界限去判斷這個值是否是我需要的。
#對每一個像元進行循環。
for i in range(rows1):
for j in range(cols1):
#b2_ds[i,j]是無效值，b3_ds[i,j]為有效值，保留b3_ds
if b2_ds[i,j] <= 0 and b3_ds[i,j] > 0:
out_ds[i, j] = b3_ds[i, j]
#b2_ds[i,j]是有效值，b3_ds[i,j]為無效值，保留b2_ds
if b3_ds[i,j] <= 0 and b2_ds[i,j] > 0:
out_ds[i,j]=b2_ds[i,j]
#b2_ds[i,j]是有效值，b3_ds[i,j]也為效值，進行均值計算
if b3_ds[i, j] > 0 and b2_ds[i, j] > 0:
out_ds[i,j]=(b2_ds[i,j]+b3_ds[i,j])/2
out_path=r'G:\0test'
outdriver = gdal.GetDriverByName('GTiff')
outds = outdriver.Create(out_path + '\\' +'AOD_mean.tif', cols1, rows1, 1,gdal.GDT_Float32)
band3:gdal.Band= outds.GetRasterBand(1)
band3.WriteArray(out_ds)
outds=None
結果如下，這裡是直接自動向上取整了。
這裡設定了背景值的忽略，是以-28672就顯示為 no data。
2.IDL中實作兩景影像的有效值合成：
這裡的測試資料是拿MCTK處理出來的MCD19A2資料，并取出了其中第二個和第三個波段進行有效值均值合成，這裡的無效值為NaN。也可以在用MCTK的過程中将填充壞值設定為一個負數。
pro mean_value
compile_opt idl2
e=envi(/headless)
raster=e.openraster('G:\000aersol\mctk_AOD.dat')
data1=raster.getdata()
;讀取栅格的行列數
nb=raster.nbands
ncols=raster.NCOLUMNS
nrows=raster.NROWS
;建立一個跟原栅格相同的數組
final_value=make_array(ncols,nrows)
for b=0,ncols-1 do begin
for c=0,nrows-1 do begin
;建立一個二維浮點型數組
data_effect=fltarr(2)
;分别讀取兩個波段的像元值
data_judge0=(data1[b,c,0] GE 0 )
data_judge1=(data1[b,c,1] GE 0 )
;将同一個像元不同兩個波段的像元值放在一個數組裡面，友善後面計算
data_effect=[data1[b,c,0],data1[b,c,1]]
;剔除掉不需要的影像值，這裡的無效值是 nan
youxiao_value=data_effect[where(data_effect ge 0)]
;進行有效值的均值合成，并将這個像元值傳給跟原栅格具有相同數組的那個數組。
final_value[b,c]=mean(youxiao_value)
endfor
endfor
out_value=final_value
uri='G:\0test\mctk_AOD_mean.dat'
final_raster=enviraster(out_value,uri=uri,spatialref=raster.spatialref)
final_raster.save
end
結果如下：
如有錯誤，歡迎指正，萬分感謝！