![a11ce56adeedd63f1256aa34049f6507.png](https://i-blog.csdnimg.cn/blog_migrate/1630fc1ac5fb6ff07617881bebf4454a.jpeg)
最近在做遥感图像时间序列云检测,用的MODIS数据。前一阵子把代码框架写完了,云检测部分暂用的产品自带的云标识。所以现在就差一个自己的云检测算法内核了。
对于一幅遥感图像,由于没有真正的标签,刚开始我用的一些非监督学习的算法(其实就是阈值法、分水岭法),阈值法是我在一幅图像上一点一点看的,可看死我了。痛定思痛,决定提取一些ROI,并把【ROI提取-格式转换-训练模型】流水线化。
1. hdf格式转换为TIFF
一个坑爹的地方是,hdf不能提ROI,因为hdf文件打开是这个样子的:
![fcf86881596ed8f22445690ab7df7996.png](https://i-blog.csdnimg.cn/blog_migrate/7cff2166c468cb3a2cb23ce0af05136b.png)
而tiff文件打开是这个样子的(我没有写入地理信息):
![082cbf00ebb0ce3cc7bbad1e55e76fbb.png](https://i-blog.csdnimg.cn/blog_migrate/86ca7cafc87917a5a7128b44f3b6d351.png)
hdf转tiff关键代码如下:
# 读入 hdf
ds = gdal.Open(hdfPath + hdf)
# 一个hdf里是有很多‘子hdf文件’的
sds = ds.GetSubDatasets()
# print(len(sds))
# for sd in sds:
# print('Name:{0}nDescription:{1}'.format(*sd))
# 读取地理和投影信息
proj = ds.GetProjection()
geoTrans = ds.GetGeoTransform()
# 从hdf中提取你需要的波段
bands = []
validBands = [13, 14, 11, 12, 15, 16, 17]
for v in validBands:
# sds[v][0]是 ‘子hdf文件’的路径, sds[v][1]是描述
band_i = gdal.Open(sds[v][0]).ReadAsArray()
bands.append(band_i)
# 写出 tiff
imgWidth, imgHeight = bands[0].shape
fileName = tiffPath+hdf[:-3]+'tiff'
# 需要先定义驱动
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(fileName, imgWidth, imgHeight, len(bands), gdal.GDT_UInt16)
# 设置地理信息
dataset.SetGeoTransform(geoTrans)
dataset.SetProjection(proj)
# 按波段写入tiff
for i in range(len(bands)):
dataset.GetRasterBand(i + 1).WriteArray(bands[i])
2. 利用ENVI提取ROI,并转换为txt格式
这一步纯用envi完成,在此不做赘述,网上有很多教程。
3. 利用pandas解析txt,将txt变成我们可用的数据集
这里主要用到了pandas库,先看txt的格式:
![3af027266473d6209ba0b652f0263bdc.png](https://i-blog.csdnimg.cn/blog_migrate/e67f2b94383352e68cad881564736835.jpeg)
思路:
1,以'n'为分隔符读入txt
2,根据第二行的number of ROIs确定数据从第几行开始
3,切割出真正的数据,再对数据进行列切割
关键代码:
txt = pd.read_csv('./w.txt', sep='n', header=None)
txt.head()
# txt.iloc[1][0]是字符串,txt.iloc[1]就是个series
numOfROI = int(txt.iloc[1][0].split()[-1])
tar = []
tarNum = []
for i in range(numOfROI):
tar.append(txt.iloc[4 + i * 4][0].split()[-1])
tarNum.append(int(txt.iloc[4 + i * 4 + 2][0].split()[-1]))
idx = 4 * (numOfROI + 1)
# 真实数据只有一列,因为在读入时是以回车为分隔符的,用.str.split(expand=True)切割成多列
data = txt.iloc[idx:, 0].str.split(expand=True)
columns = txt.iloc[idx-1, 0].split()[1:]
data.columns = columns
# data.rename(columns=dict(zip(data.columns, columns)), inplace=True)
# 删除掉不需要的列
data.drop(['ID', 'X', 'Y'], axis=1, inplace=True)
data.reset_index(drop=True,inplace=True)
# 增加标签列
y = np.zeros(sum(tarNum), np.int8)
y[tarNum[0]:] = 1
data['class'] = y
# 注意不要保存索引
data.to_csv('./w.csv', index=False)
4. 利用sklearn库训练多个模型
有了数据集,接下来的事就很简单了。将多个模型存入到一个字典中,直接训练多个模型比较结果。
#%%
import pandas as pd
from sklearn import svm
from sklearn.model_selection import train_test_split
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import RidgeCV
from sklearn.metrics import accuracy_score
#%%
df = pd.read_csv('./w.csv')
col = df.columns
X_train, X_test, y_train, y_test = train_test_split(df[col[:-1]], df[col[-1]],
test_size=0.33, random_state=42, stratify=df[col[-1]])
ESTIMATORS = {
"Extra trees": ExtraTreesRegressor(n_estimators=10, max_features=4,
random_state=0),
"K-nn": KNeighborsRegressor(),
"Linear regression": LinearRegression(),
"Ridge": RidgeCV(),
"SVM": svm.SVC(kernel='poly')
}
y_test_predict = dict()
for name, estimator in ESTIMATORS.items():
estimator.fit(X_train, y_train)
y_test_predict[name] = estimator.predict(X_test)
print('{0:20}: {1}'.format(name, accuracy_score(y_test, y_test_predict[name]>0.5)))
发现一个有趣的结果:
![0076153c7fa4a12d123c1b571dfdbc6b.png](https://i-blog.csdnimg.cn/blog_migrate/5acf64970e40dcf445212b051325f568.png)
![1ec003f04bf656ee44338531bf210799.png](https://i-blog.csdnimg.cn/blog_migrate/df8379629e96b40026b63f939232510e.png)
第一个svm用的是rbf核,效果很不好,在我发现ET模型有很好的效果后,将SVM改为多项式核,果然取得了很好的效果。树模型其实就是阈值法,对数据进行加减法,直觉上感觉和多项式很像。