囚徒4.0_11_基于python的风云云检测算法

#囚徒4.0_11_基于python的风云算法
#关于昨天数据不同的问题：是因为IDL和Python的逻辑不同而导致的，数据读取没问题，我表示错了。
#换语言好麻烦，现在都不知道什么语法对应什么语言了，一团糟。
#从上午十点写到现在，测试的时候发现python他的读取逻辑和IDL不一样，他的循环也不一样，我真的栓Q了。
from itertools import count
import h5py
import numpy as np
from osgeo import gdal
from pyhdf.SD import SD, SDC
#风云4km地理数据读取 正确
def fy_4km_geo_read(file):
with h5py.File(file, 'r') as f:
# 可以查看文件中包含的所有数据集名称
#print(list(f.keys()))
SZA4km=f['/Navigation☺MSunZenith'][:]
dims4km=SZA4km.shape
#print('里面: SZA(1400,820)的值:{}'.format(SZA4km[1400,820]))
return SZA4km,dims4km
#风云4km波段数据读取 正确
def fy_4km_fid_read(file,SZA,dims4km):
#***1 数据准备
DEGRAD=np.arccos(-1.)/180.0 #将180度转换成弧度，得到弧度与角度之间的转换系数，存储在变量DEGRAD中。
Zero_Eps=0.000001 #定义一个很小的常数，用来比较浮点数是否为0。
Rel_Equality_Eps=0.000001 #定义一个很小的常数，用来判断两个浮点数是否相等。
SolZen_Threshold=84.0 #定义太阳天顶角的阈值，用来判断哪些数据点的反射率是有效的。
FV_L1B = -999.0 #定义一个无效值，用来标记无法计算出TOA反射率的数据点。
cossza=np.cos(DEGRAD*SZA) #计算太阳天顶角的余弦值，并将结果存储在变量cossza中。
cos_ind = np.where(abs(cossza) >= Zero_Eps, 1, 0) #根据预定义的常数，判断太阳天顶角余弦值的绝对值是否大于等于一个很小的正数Zero_Eps，并将结果存储在变量cos_ind中。
x_cossza = (1.0/cossza)*cos_ind #计算一个系数x_cossza，用来修正反射率数据。
geo_valid =np.where((SZA > 0)&(SZA < SolZen_Threshold)& (cos_ind == 1), 1, 0)
#根据预定义的条件，生成有效的反射率数据。这个表达式可以看作是两部分的逻辑与操作，第一部分是太阳天顶角在0到阈值之间，第二部分是太阳天顶角的余弦值在一个很小的正数以上。只有同时满足这两个条件，才能判定该点的反射率数据有效。
#波段数据读取
with h5py.File(file, 'r') as f:
#波段读取
nb1=f['/Data☺MChannel01'][:]
nb2=f['/Data☺MChannel02'][:]
nb3=f['/Data☺MChannel03'][:]
nb4=f['/Data☺MChannel04'][:]
nb5=f['/Data☺MChannel05'][:]
nb6=f['/Data☺MChannel06'][:]
nb7=f['/Data☺MChannel07'][:]
nb8=f['/Data☺MChannel08'][:]
nb9=f['/Data☺MChannel09'][:]
nb10=f['/Data☺MChannel10'][:]
nb11=f['/Data☺MChannel11'][:]
nb12=f['/Data☺MChannel12'][:]
nb13=f['/Data☺MChannel13'][:]
nb14=f['/Data☺MChannel14'][:]
nb15=f['/Data☺MChannel15'][:]
#定标表
cal_nb1=f['/Calibration/CALChannel01'][:]
cal_nb2=f['/Calibration/CALChannel02'][:]
cal_nb3=f['/Calibration/CALChannel03'][:]
cal_nb4=f['/Calibration/CALChannel04'][:]
cal_nb5=f['/Calibration/CALChannel05'][:]
cal_nb6=f['/Calibration/CALChannel06'][:]
cal_nb7=f['/Calibration/CALChannel07'][:]
cal_nb8=f['/Calibration/CALChannel08'][:]
cal_nb9=f['/Calibration/CALChannel09'][:]
cal_nb10=f['/Calibration/CALChannel10'][:]
cal_nb11=f['/Calibration/CALChannel11'][:]
cal_nb12=f['/Calibration/CALChannel12'][:]
cal_nb13=f['/Calibration/CALChannel13'][:]
cal_nb14=f['/Calibration/CALChannel14'][:]
cal_nb15=f['/Calibration/CALChannel15'][:]
cal_nb1=np.append(cal_nb1,-999.0) #后面添加一个标记值，用于记录各波段的无效值。
cal_nb2=np.append(cal_nb2,-999.0)
cal_nb3=np.append(cal_nb3,-999.0)
cal_nb4=np.append(cal_nb4,-999.0)
cal_nb5=np.append(cal_nb5,-999.0)
cal_nb6=np.append(cal_nb6,-999.0)
cal_nb7=np.append(cal_nb7,-999.0)
cal_nb8=np.append(cal_nb8,-999.0)
cal_nb9=np.append(cal_nb9,-999.0)
cal_nb10=np.append(cal_nb10,-999.0)
cal_nb11=np.append(cal_nb11,-999.0)
cal_nb12=np.append(cal_nb12,-999.0)
cal_nb13=np.append(cal_nb13,-999.0)
cal_nb14=np.append(cal_nb14,-999.0)
cal_nb15=np.append(cal_nb15,-999.0)
#***3 计算反射率
ref_b1=cal_nb1[np.minimum(nb1, cal_nb1.size-1)]
ref_valid=np.where((ref_b1 > 0.0)&(ref_b1 <1.0)&(geo_valid == 1),1,0)
ref_b1=(ref_b1*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
ref_b2=cal_nb2[np.minimum(nb2, cal_nb2.size-1)] #根据读取到的数据和波段编号，计算出指定波段的TOA反射率，将其存储在变量ref_b1中
ref_valid=np.where((ref_b2 > 0.0)&(ref_b2 <1.0)&(geo_valid == 1),1,0)
ref_b2=(ref_b2*ref_valid)*x_cossza + FV_L1B*(1-ref_valid) #对计算得到的TOA反射率进行修正和校正，得到最终的输出结果。这个表达式可以看作是两部分的加权和，第一部分是计算得到的TOA反射率乘以一个系数x_cossza，第二部分是一个无效值FV_L1B乘以另一个系数1-temporary(ref_valid)。只有在TOA反射率有效时，才需要对其进行修正和校正。
ref_b3=cal_nb3[np.minimum(nb3, cal_nb3.size-1)]
ref_valid=np.where((ref_b3 > 0.0)&(ref_b3 <1.0)&(geo_valid == 1),1,0)
ref_b3=(ref_b3*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
ref_b4=cal_nb4[np.minimum(nb4, cal_nb4.size-1)]
ref_valid=np.where((ref_b4 > 0.0)&(ref_b4 <1.0)&(geo_valid == 1),1,0)
ref_b4=(ref_b4*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
ref_b5=cal_nb5[np.minimum(nb5, cal_nb5.size-1)]
ref_valid=np.where((ref_b5 > 0.0)&(ref_b5 <1.0)&(geo_valid == 1),1,0)
ref_b5=(ref_b5*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
ref_b6=cal_nb6[np.minimum(nb6, cal_nb6.size-1)]
ref_valid=np.where((ref_b6 > 0.0)&(ref_b6 <1.0)&(geo_valid == 1),1,0)
ref_b6=(ref_b6*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
#亮温
tmp_7=cal_nb7[np.minimum(nb7, cal_nb7.size-1)]
ref_valid=np.where((tmp_7 != -999),1,0)
tmp_7=(tmp_7*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_8=cal_nb8[np.minimum(nb8, cal_nb8.size-1)]
ref_valid=np.where((tmp_8 != -999),1,0)
tmp_8=(tmp_8*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_9=cal_nb9[np.minimum(nb9, cal_nb9.size-1)]
ref_valid=np.where((tmp_9 != -999),1,0)
tmp_9=(tmp_9*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_10=cal_nb10[np.minimum(nb10, cal_nb10.size-1)]
ref_valid=np.where((tmp_10 != -999),1,0)
tmp_10=(tmp_10*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_11=cal_nb11[np.minimum(nb11, cal_nb11.size-1)]
ref_valid=np.where((tmp_11 != -999),1,0)
tmp_11=(tmp_11*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_12=cal_nb12[np.minimum(nb12, cal_nb12.size-1)]
ref_valid=np.where((tmp_12 != -999),1,0)
tmp_12=(tmp_12*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_13=cal_nb13[np.minimum(nb13, cal_nb13.size-1)]
ref_valid=np.where((tmp_13 != -999),1,0)
tmp_13=(tmp_13*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_14=cal_nb14[np.minimum(nb14, cal_nb14.size-1)]
ref_valid=np.where((tmp_14 != -999),1,0)
tmp_14=(tmp_14*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
tmp_15=cal_nb15[np.minimum(nb15, cal_nb15.size-1)]
ref_valid=np.where((tmp_15 != -999),1,0)
tmp_15=(tmp_15*ref_valid)*x_cossza + FV_L1B*(1-ref_valid)
#fy_4km_fid_read end
return ref_b2,ref_b3,ref_b4,ref_b5,tmp_7,tmp_8,tmp_9,tmp_10,tmp_13,tmp_15
#掩膜数据读取。注意路径不能为中文，pyhdf库不支持with读取文件，本文件为hdf4格式 正确
def lan_sea_mask_read(file):
f=SD(file, SDC.READ)
#print(f.datasets().keys())
HGHT = f.select('Height')[:]
LSM = f.select('Land/SeaMask')[:]
#海陆掩膜优化
mask_lan=np.where((LSM == 1)|(LSM == 2)|(LSM == 4),1,0)
ns=mask_lan.shape
# 0：浅海洋（靠近海岸线距离小于5千米或者深度小于50米的海洋区域）。
# 1：陆地（除了其他所有类别以外的土地区域）。
# 2：海洋海岸线和湖泊沿岸线。
# 3：浅内陆水体（位于海岸线距离小于5千米或者深度小于50米的内陆水体）。
# 4：短暂性水体（间歇性的水体）。
# 5：深内陆水体（位于海岸线距离大于5千米并且深度大于50米的内陆水体）。
# 6：温带或者大陆型海洋（位于距离海岸线大于5千米、深度大于50米且小于500米的海洋区域）。
# 7：深海洋（海洋深度大于500米）
#遍历全部，当5*5窗口内的陆地像元数大于窗口像元的0.6则判定此像元为陆地像元
W=2
for i in range(W,ns[0]-W):
for j in range(W,ns[1]-W):
if mask_lan[i,j]==0:
count=0
#窗口
for m in range(i-W,i+W+1):
for n in range(j-W,j+W+1):
if mask_lan[m,n]== 1:
count+=1
if count >15:
mask_lan[i,j]=1
#lan_sea_mask_read end
return mask_lan,HGHT
#可信矩阵编写 正确
def T_mat(band,cloud,cloud_f,iden):
sz=band.shape
temp=np.where(band!= -999,1,0)
trust=np.zeros((sz[0],sz[1]),dtype=float)
if iden=='gt':
for i in range(0,sz[0]):
for j in range(0,sz[1]):
if (band[i,j]> cloud)&(temp[i,j]==1):
trust[i,j]=0
elif (band[i,j]< cloud_f)&(temp[i,j]==1):
trust[i,j]=1
else:
trust[i,j]=(cloud-band[i,j])/(cloud-cloud_f)*1.0
else:
for i in range(0,sz[0]):
for j in range(0,sz[1]):
if (band[i,j]< cloud)&(temp[i,j]==1):
trust[i,j]=0
elif (band[i,j]> cloud_f)&(temp[i,j]==1):
trust[i,j]=1
else:
trust[i,j]=(band[i,j]-cloud)/(cloud_f-cloud)*1.0
#T_mat end
return trust
#可信度插值
def inter(trust_mat,x,y,ratio,scale):
count=0
tot=0
if ratio>=0.5:
for i in range(x-1,x+2):
for j in range(y-1,y+2):
if trust_mat[i,j]>0.5:
tot+=trust_mat[i,j]*ratio
count+=1
else:
tot+=trust_mat[i,j]*(1-ratio)
value=tot/(ratio*count+(1-ratio)*(9-count))
value*=scale
else:
for i in range(x-1,x+2):
for j in range(y-1,y+2):
if trust_mat[i,j]<0.5:
tot+=trust_mat[i,j]*ratio
count+=1
else:
tot+=trust_mat[i,j]*(1-ratio)
value=tot/(ratio*count+(1-ratio)*(9-count))
value*=scale
#inter end
return value
#海陆边线
def math_lan(mask_lan):
sz=mask_lan.shape
lan_line=np.zeros((sz[0],sz[1]))
W=1
for i in range(1,sz[0]-1):
for j in range(1,sz[1]-1):
mask_0=0
mask_1=1
for m in range(i-W,i+W+1):
for n in range(j-W,j+W+1):
if mask_lan [m,n]==0:
mask_0=0
else:
mask_1=1
if (mask_0 == 0)&(mask_1== 1):
lan_line[i,j]=1
#math_lan end
return lan_line