c#写图像tif gdal_遥感图像分类智能化

最近在做遥感图像时间序列云检测，用的MODIS数据。前一阵子把代码框架写完了，云检测部分暂用的产品自带的云标识。所以现在就差一个自己的云检测算法内核了。

对于一幅遥感图像，由于没有真正的标签，刚开始我用的一些非监督学习的算法（其实就是阈值法、分水岭法），阈值法是我在一幅图像上一点一点看的，可看死我了。痛定思痛，决定提取一些ROI，并把【ROI提取-格式转换-训练模型】流水线化。

1. hdf格式转换为TIFF

一个坑爹的地方是，hdf不能提ROI，因为hdf文件打开是这个样子的：

而tiff文件打开是这个样子的（我没有写入地理信息）：

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的格式：

思路：

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)))

发现一个有趣的结果：

第一个svm用的是rbf核，效果很不好，在我发现ET模型有很好的效果后，将SVM改为多项式核，果然取得了很好的效果。树模型其实就是阈值法，对数据进行加减法，直觉上感觉和多项式很像。