SVM检测影像的云

首先利用arcgis对landsat影像打点云样本和非云样本点图层，转为csv，用作机器学习检测样本

import numpy as np

import pandas as pd

from osgeo import gdal, gdal_array

# Tell GDAL to throw Python exceptions, and register all drivers

gdal.UseExceptions()

#gdal.AllRegister()#不知道为啥总是报错

# Read in raster image

img_ds = gdal.Open('D:/idltemp/yw/2014/cloubtest.tif', gdal.GA_ReadOnly)

img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount),

gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType))

for b in range(img.shape[2]):

img[:, :, b] = img_ds.GetRasterBand(b + 1).ReadAsArray()

new_shape = (img.shape[0] * img.shape[1], img.shape[2])

print (img.shape)

print (new_shape)

X = img[:, :, :3].reshape(new_shape)

from sklearn import svm

datacloudtest=pd.read_csv(r'D:\idltemp\yw\2014\1.csv')

datacloudtest=np.array(datacloudtest)

cloud,band1,band2,band3=datacloudtest[:,1],datacloudtest[:,2],datacloudtest[:,3],datacloudtest[:,4]

cloudtest=np.column_stack((cloud,band1,band2,band3))

ytest=cloudtest[:,0].reshape(-1,1).astype('int')

xtest=cloudtest[:,1:]

cloudy_train = ytest[:int(len(ytest)*0.9)]#设置训练样本和测试样本

cloudy_test = ytest[int(len(ytest)*0.9):]

#y_test = y2[int(len(y2)*0.9):]

cloudX_train = xtest[:int(len(ytest)*0.9)]

cloudX_test = xtest[int(len(ytest)*0.9):]

#cloudgbm = lgb.LGBMRegressor(boosting_type='gbdt', num_leaves=31, max_depth=-1, learning_rate=0.1, n_estimators=100, subsample_for_bin=200000, objective=None, class_weight=None, min_split_gain=0.0, min_child_weight=0.001, min_child_samples=20, subsample=1.0, subsample_freq=0, colsample_bytree=1.0, reg_alpha=0.0, reg_lambda=0.0, random_state=None, n_jobs=-1, silent=True, importance_type='split',)#参数设置

cloudgbm =svm.SVC(gamma='scale')

cloudgbm.fit(cloudX_train, cloudy_train)

cloudy_pred = cloudgbm.predict(X)

cloudy_pred = cloudy_pred.reshape(img[:, :, 0].shape)

import matplotlib.pyplot as plt

print (cloudy_pred.shape)

plt.figure(figsize=(20,20))

plt.imshow(cloudy_pred, cmap="hsv")

plt.show()

from osgeo import gdal, gdal_array

import gdal

## write out to tiff

### single band raster.

ds = gdal.Open('D:/idltemp/yw/2014/cloubtest.tif')

band = ds.GetRasterBand(2)

arr = band.ReadAsArray()

x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()

[cols, rows] = arr.shape

format = "GTiff"

driver = gdal.GetDriverByName(format)

outDataRaster = driver.Create("D:/idltemp/yw/2014/cloudpre1.tif", rows, cols, 1, gdal.GDT_Int16)#如果已经存在文件，也会runtimeerror

outDataRaster.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input

outDataRaster.SetProjection(ds.GetProjection())##sets same projection as input

outDataRaster.GetRasterBand(1).WriteArray(cloudy_pred)

outDataRaster.FlushCache() ## remove from memory...important

del outDataRaster ## delete the data (not the actual geotiff)

结果

比较