Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model，ESTARFM

前言

面对华东窗口这类问题，必须要有的一些基础知识

• xy2index
• define search extents

Paper Reading

OK

Code Interpreting

IDL

Unread

Python

Continuously updating……

# utils.py
# 使用numpy和gdal库写了2个函数，分别实现影像的读和写。
# ------------------------------------------------------------

# parameters_estarfm.yaml

w: 99  # set the half window size, if 25, the window size is 25*2+1=51 fine pixels
num_class: 6.0  # set the estimated number of classes, please set a larger value if blending images with very few bands
DN_min: 0.0  # set the range of DN value of the image,If byte, 0 and 255
DN_max: 10000.0
background: -9999   # set the value of background pixels. 0 means that pixels will be considered as background if one of its bands= 0
patch_long: 500   # set the size of each block,if process whole ETM scene, set 500-1000

Digital Number, DN :https://baike.baidu.com/item/DN值/8633825?fr=aladdin

# ESTARFM.py

# import gdal 改为 from osgeo import gdal
# np.int 改为 np.int_

1 读取参数
2 设置输出路径
3 读取五幅影像
3.1 将影像切块（向上取整获得块的行列数）
3.2 切块（block）后的影像保存到临时文件夹

------------------------------ 下面开始对每个块block进行操作 --------------------------------------------------------------------------------------------------
3.3 将切好块的影像依次读取（第1个block的5个影像f1c1f2c2c0、第二个block的五个影像…… 并为预测影像分配空间
3.4 相似像元阈值

# compute the threshold of similar pixel seeking
similar_th = np.zeros((2, nb)).astype(float)
for iband in range(0, nb):
    similar_th[0, iband] = np.std(fine1[iband, :, :] * 2.0 / num_class)  # 第一幅精细影像的阈值
    similar_th[1, iband] = np.std(fine2[iband, :, :] * 2.0 / num_class)  # 第二幅精细影像的阈值

3.5 利用idlwrap中的索引计算欧氏距离后再转换为geographic distance

# compute the distance of each pixel in the window with the target pixel (integrate window)
D_temp1 = w - np.tile((idlwrap.indgen(w*2+1)), (int(w*2+1), 1))  # (1 2 3 …… 2*w+1)|(w*2+1,)*w*2+1->(w*2+1, 2*2+1)
d1 = np.power(D_temp1, 2)  # 乘方
D_temp2 = w - np.tile(idlwrap.indgen(1, w*2+1), (1, int(w*2+1)))
d2 = np.power(D_temp2, 2)
D_D_all = 1.0 + np.sqrt(d1 + d2) / float(w)
D_D_all = D_D_all.flatten()

3.6 筛选背景值
有效像元标记为1，
-----------------------------------------------------------下边是在每一个块中使用滑动窗口进行操作-------------------------------------------------------------
nl，ns是每个block的行列号，
3.7 确定滑动窗口的范围
3.8 搜索相似像元的值
文章中提到，在同一像元内的像元有着相同的转换系数，所以将相似像元选中后展品成向量便于后续计算
3.9 计算相关系数
如果筛选出来的，符合条件的相似像元数量大于5
3.10

测试影像运行时间

---

这个代码很多地方可以改进，我将会对内部的一些函数做进一步封装，利用GUP做加速……