IDL&python学习——实现两景影像中有效值的均值合成

最新推荐文章于 2024-12-05 15:11:21 发布

置顶三十二号星期八

最新推荐文章于 2024-12-05 15:11:21 发布

阅读量2.3k

点赞数 8

分类专栏： python 有效值均值合成 IDL 文章标签： python 其他

本文链接：https://blog.youkuaiyun.com/weixin_42176976/article/details/108143435

版权

IDL 同时被 3 个专栏收录

10 篇文章

订阅专栏

python

6 篇文章

订阅专栏

有效值均值合成

1 篇文章

订阅专栏

最近在处理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

结果如下：