随机森林

库的知识补充

1、pd.isna（大部分数据处理的场合推荐使用）

判断DateFrame格式的数据是否出现空值，如果我们的空值既会出现np.nan，也会出现math.nan，甚至还会出现None，或者要判断的数据值既可能是数值型也可能是字符串（符合大部分场景的实际情况），墙裂推荐使用pd.isna，例如下面代码：

pd.isna(np.nan)
Out[29]: True
pd.isna(math.nan)
Out[30]: True
pd.isna(None)
Out[31]: True
pd.isna('a')
Out[32]: False
pd.isna(10)
Out[33]: False

当测试或实际应用需要设置空值缺失值时可以用None 、np.nan 、pd.NaT
处理空值异常值有两种办法**，一种是将空值行/列删除，一种是将替代空值**
如何对空值计数
a=df.isnull()
b=a[a==True]
b.count()#用来计算nan数量

删除的方法
将nan的行全部删除
df.dropna()
print(‘dropna’,df.dropna())#将带nan的行删除 axis=1删除列
或者是 将某列行不为空值或者某几列行均不为空值的行另存为
data = data[data[‘SUM_YR_1’].notnull()*data[‘SUM_YR_2’].notnull()]
填充
df.fillna(10) 自动按10填充
df.Age.fillna(df.Age.mean(),inplace = True) 用平均值来填充

print(‘ffill’,df.fillna(method=‘ffill’))#向上填充
print(‘dfill’,df.fillna(method=‘bfill’))#向下填充
print(‘replace’,df.replace({np.nan:‘aa’}))#对应的填充替换

2、np.array格式的数据处理

print(np.isnan(data).any())，判断数据是否有空值
data[np.isnan(data)] = 0，将控制赋值为0
Numpy数组(ndarray)中含有缺失值(nan)行和列的删除方法

print(np.isnan(a).any(axis=1))
print(~np.isnan(a).any(axis=1))
print(a[~np.isnan(a).any(axis=1), :]
print(a[~np.isnan(a).any(axis=1)])

3、pickle模型保存与加载

pickle.dumps()将对象obj对象序列化并返回一个byte对象
pickle.loads(),从字节对象中读取被封装的对象

4、python字典怎么根据值返回键

随机森林的函数

class sklearn.ensemble.RandomForestClassifier(n_estimators=10, criterion='gini', max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=’auto’, max_leaf_nodes=None, min_impurity_decrease=0.0, min_impurity_split=None, bootstrap=True, oob_score=False, n_jobs=1, random_state=None, verbose=0, warm_start=False, class_weight=None)

参数：
sklearn随机森林分类类RandomForestClassifier
n_estimators : integer, optional (default=10) 整数，可选择(默认值为100)决策数的数目
criterion : string, optional (default=”gini”) 字符串，可选择(默认值为“gini”)。特征选择原则
max_features : int, float, string or None, optional (default=”auto”) 整数，浮点数，字符串或者无值，可选的（默认值为"auto"）

代码

1、遥感影像的读取与处理

1、 读取数据到TXT

#读取tif数据集
def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName+"文件无法打开")
    return dataset
dataset = readTif(Landset_Path)  #首先将栅格图像读取
'''
获取栅格影像的基本信息
'''
@author: User
将遥感影像像素值和标签值读取到txt文件中
"""

import gdal
import os
import random

#读取tif数据集
def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName+"文件无法打开")
    return dataset


Landset_Path = r"E:\LZ\ROI\growth.tif"  #遥感图像
LabelPath = r"E:\LZ\ROI\train_roi\train_set_4800_clip.tif"  #标签图片
txt_Path=r"E:\LZ\Data\data1.txt"  #存放数的文件

##########################################################  1、读取图像数据
dataset = readTif(Landset_Path)  #首先将栅格图像读取
'''
获取栅格影像的基本信息
'''
Tif_width = dataset.RasterXSize #栅格矩阵的列数
Tif_height = dataset.RasterYSize #栅格矩阵的行数
Tif_bands = dataset.RasterCount #波段数
Tif_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息
Tif_proj = dataset.GetProjection()#获取投影信息
Landset_data = dataset.ReadAsArray(0,0,Tif_width,Tif_height)#读取影像值

dataset = readTif(LabelPath)#2、读取标注信息
Label_data = dataset.ReadAsArray(0,0,Tif_width,Tif_height)#获取标注信息的数值


# 写之前，先检验文件是否存在，存在就删掉
if os.path.exists(txt_Path):
    os.remove(txt_Path)
# 以写的方式打开文件，如果文件不存在，就会自动创建
file_write_obj = open(txt_Path, 'w') #open的方式打开文件，需要后面关闭文件。而with open不需要

####################################################首先收集植被类别样本，
####################################################遍历所有像素值，
####################################################为植被的像元全部收集。
count = 0
for i in range(Label_data.shape[0]):#行
    for j in range(Label_data.shape[1]):#列
        #  我设置的植被类别在标签图中像元值为1
        if(Label_data[i][j] == 1):
            var = ""
            for k in range(Landset_data.shape[0]):#波段
                var = var + str(Landset_data[k][i][j])+","
            var = var + "surgarcane"
            file_write_obj.writelines(var)
            file_write_obj.write('\n')
            count = count + 1


count = 0
for i in range(Label_data.shape[0]):
    for j in range(Label_data.shape[1]):
        #  我设置的植被类别在标签图中像元值为1
        if(Label_data[i][j] == 2):
            var = ""
            for k in range(Landset_data.shape[0]):
                var = var + str(Landset_data[k][i][j])+","
            var = var + "forest"
            file_write_obj.writelines(var)
            file_write_obj.write('\n')
            count = count + 1   
            
            

count = 0
for i in range(Label_data.shape[0]):
    for j in range(Label_data.shape[1]):
        #  我设置的植被类别在标签图中像元值为1
        if(Label_data[i][j] == 3):
            var = ""
            for k in range(Landset_data.shape[0]):
                var = var + str(Landset_data[k][i][j])+","
            var = var + "water"
            file_write_obj.writelines(var)
            file_write_obj.write('\n')
            count = count + 1   


count = 0
for i in range(Label_data.shape[0]):
    for j in range(Label_data.shape[1]):
        #  我设置的植被类别在标签图中像元值为1
        if(Label_data[i][j] == 4):
            var = ""
            for k in range(Landset_data.shape[0]):
                var = var + str(Landset_data[k][i][j])+","
            var = var + "buildings"
            file_write_obj.writelines(var)
            file_write_obj.write('\n')
            count = count + 1  


count = 0
for i in range(Label_data.shape[0]):
    for j in range(Label_data.shape[1]):
        #  我设置的植被类别在标签图中像元值为1
        if(Label_data[i][j] == 5):
            var = ""
            for k in range(Landset_data.shape[0]):
                var = var + str(Landset_data[k][i][j])+","
            var = var + "others"
            file_write_obj.writelines(var)
            file_write_obj.write('\n')
            count = count + 1  
           
####################################################其次收集非植被类别样本，
####################################################因为非植被样本比植被样本多很多，
####################################################所以采用在所有非植被类别中随机选择非植被样本，
####################################################数量与植被样本数量保持一致。
#Threshold = count
#count = 0
#for i in range(10000000000):
#    X_random = random.randint(0,Label_data.shape[0]-1)
#    Y_random = random.randint(0,Label_data.shape[1]-1)
#    #  我设置的非植被类别在标签图中像元值为0
#    if(Label_data[X_random][Y_random] == 0):
#        var = ""
#        for k in range(Landset_data.shape[0]):
#            var = var + str(Landset_data[k][X_random][Y_random])+","
#        var = var + "Unclassified"
#        file_write_obj.writelines(var)
#        file_write_obj.write('\n')
#        count = count + 1
#    if(count == Threshold):
#        break

file_write_obj.close()

2、训练模型

from sklearn.ensemble import RandomForestClassifier
import numpy as np
from sklearn import model_selection
import pickle 
from sklearn.model_selection import KFold
from sklearn.metrics import recall_score #二分类问题的召回率
from sklearn.metrics import cohen_kappa_score #多分类问题的kappa系数
import pandas as pd
##  定义字典，便于来解析样本数据集txt
#def Iris_label(s):
#    it={b'Vegetation':0, b'Non-Vegetation':1}
#    return it[s]

def Iris_label(s):
    it={b'surgarcane':1,b'forest':2,b'water':3,b'buildings':4,b'others':5}
    return it[s]

path=r"E:\LZ\Data\data.txt"
SavePath = r"E:\LZ\Data\model.pickle"


#  1.读取数据集
#data=np.loadtxt(path, dtype=float, delimiter=',', converters={7:Iris_label} )
data=np.loadtxt(path, dtype=float, delimiter=',', converters={4:Iris_label} )
#  converters={7:Iris_label}中“7”指的是第8列：将第8列的str转化为label(number)
#2.划分数据与标签
x,y=np.split(data,indices_or_sections=(4,),axis=1) #x为数据，y为标签,axis纵向切分
print(y)
#x=x[:,0:4] #选取前7个波段作为特征
train_data,test_data,train_label,test_label = model_selection.train_test_split(x,y, random_state=1, train_size=0.9,test_size=0.1)
print(train_label)
print(test_label)
#3、对空值数据进行赋值处理
train_data[np.isnan(train_data)] = 0 #数据中有空值，将其赋值为0
test_data[np.isnan(test_data)] = 0
#数据进行预处理，根据情况选定
mean_image = np.mean(train_data, axis=0)
train_data-= mean_image
test_data -= mean_image

#  4.用100个树来创建随机森林模型，训练随机森林
classifier = RandomForestClassifier(n_estimators=100, 
                               bootstrap = True,
                               max_features = 'sqrt')
classifier.fit(train_data, train_label.ravel())#ravel函数拉伸到一维
#  5.计算随机森林的准确率
print("训练集：",classifier.score(train_data,train_label))
print("测试集：",classifier.score(test_data,test_label))
#
#  6.保存模型
#以二进制的方式打开文件：
file = open(SavePath, "wb")
#将模型写入文件：
pickle.dump(classifier, file)#pickle.dumps()将对象obj对象序列化并返回一个byte对象
#最后关闭文件：
file.close()

模型参数调优采用K折交叉验证法
注意：在sklearn.metrics 模块中 recall_score 和 precision_score函数不支持计算多分类
多标签分类的评价指标
多标签学习系统中的性能评估与经典的单标签学习问题中不同，在单标签问题中使用的经典度量标准包括准确率，Precision，Recall 和 F-measure.在多标签学习中，评估会更加复杂，对于一个测试集。
评价准则1：Kappa系数
Kappa系数是基于混淆矩阵的计算得到的模型评价参数。计算公式如下

系数的值在-1到1之间，系数小于0的话实际上就相当于随机了。python实现为:

系数的值在-1到1之间，系数小于0的话实际上就相当于随机了。
python实现为:

from sklearn.metrics import cohen_kappa_score
kappa = cohen_kappa_score(y_true,y_pred，label=None) #(label除非是你想计算其中的分类子集的kappa系数，否则不需要设置)

2.海明距离
海明距离也适用于多分类的问题，简单来说就是衡量预测标签与真实标签之间的距离，取值在0~1之间。距离为0说明预测结果与真实结果完全相同，距离为1就说明模型与我们想要的结果完全就是背道而驰。python实例。

from sklearn.metrics import hamming_loss
ham_distance = hamming_loss(y_true,y_pred)

3.杰卡德相似系数
它与海明距离的不同之处在于分母。当预测结果与实际情况完全相符时，系数为1；当预测结果与实际情况完全不符时，系数为0；当预测结果是实际情况的真子集或真超集时，距离介于0到1之间。
我们可以通过对所有样本的预测情况求平均得到算法在测试集上的总体表现情况。
from sklearn.metrics import jaccard_similarity_score
jaccrd_score = jaccrd_similarity_score(y_true,y_pred,normalize = default)
#normalize默认为true，这是计算的是多个类别的相似系数的平均值，normalize = false时分别计算各个类别的相似系数

4.铰链损失
铰链损失（Hinge loss）一般用来使“边缘最大化”（maximal margin）。损失取值在0~1之间，当取值为0，表示多分类模型分类完全准确，取值为1表明完全不起作用。
from sklearn.metrics import hinge_loss
hinger = hinger_loss(y_true,y_pred)

#K折交叉验证，对参数进行调优
#为进行交叉验证转化数据格式
train_data=pd.DataFrame(train_data)
train_label=pd.DataFrame(train_label)
def printing_Kfold_scores(x_train_data,y_train_data):
    #k折交叉验证
    fold = KFold(n_splits=5,shuffle=False)
    #不同的C参数
#    c_param_range = [0.01,0.1,1,10,100]
    
#    results_table = pd.DataFrame(index = range(len(c_param_range),2), columns = ['C_parameter','Mean recall score'])
#    results_table['C_parameter'] = c_param_range
    
    n_estimators_param_range = [50,70,90,100,110,120]  #随机森林的树的数目
    max_depth_param_range=[2,3,4,5]#每棵树的最大深度
    
#    results_table = pd.DataFrame(index = range(len(n_estimators_param_range),2), columns = ['n_estimators_parameter','Mean recall score'])
#    results_table['n_estimators_parameter'] = c_param_range
    #k折操作将会给出两个列表：train_indices = indices[0], test_indices = indices[1]
    results = {}
    best_mean_recall_score = 0
    best_net = None
    print("-----------starting trainging-------------")
    for ne in  n_estimators_param_range:
        for md in max_depth_param_range:
            recall_accs = []
            train_scores=[]
            test_scores=[]
            for it,ids in enumerate(fold.split(x_train_data),start=1):
                
                classifier = RandomForestClassifier(n_estimators=ne, 
                                                    max_depth=md,
                                                    bootstrap = True,
                                                    max_features = 'sqrt')
                classifier.fit(x_train_data.iloc[ids[0]],y_train_data.iloc[ids[0]].values.ravel())
                y_pred_undersample = classifier.predict(x_train_data.loc[ids[1],:].values)
#                print(y_train_data.loc[ids[1],:].values)
#                recall_acc = recall_score(y_train_data.loc[ids[1],:].values,y_pred_undersample.ravel())#召回率只能处理二分类问题
#                recall_acc = np.mean(y_train_data.loc[ids[1],:].values==y_pred_undersample.ravel())#此次先采用准确率
                kappa = cohen_kappa_score(y_train_data.loc[ids[1],:].values,y_pred_undersample.ravel())
                recall_accs .append(kappa)
                print('Iteration ', it,': recall score = ', kappa)
                train_score=classifier.score(x_train_data.loc[ids[0]],y_train_data.loc[ids[0]].values)
                train_scores.append(train_score)
#                print("训练集：",train_scores)
                test_score=classifier.score(x_train_data.loc[ids[1],:].values,y_train_data.loc[ids[1]].values)
                test_scores.append(test_score)
#                print("测试集：",train_scores)
            mean_recall_score =np.mean(recall_accs)#计算的平均召回率
#            mean_train_scores=np.mean(recall_accs)
#            mean_test_scores=np.mean(test_scores)
             
            print('Mean recall score ', mean_recall_score)
            #求出最大的召回率，何其所对应的分类器，将参数结果保存到results中
            if mean_recall_score> best_mean_recall_score :
                best_mean_recall_score =mean_recall_score
                best_net=classifier
                
            results[(ne,md)]= mean_recall_score
    for ne, md in sorted(results):
        val_acc = results[(ne, md )]
        print ('ne %d md %d  val accuracy: %f' % (ne, md,val_acc))
    best_params=list(results.keys())[list(results.values()).index(best_mean_recall_score)]
    print( 'best validation accuracy achieved during cross-validation: %f' %  best_mean_recall_score)
    print( 'best params achieved during cross-validation: %f' ,best_params)
    return best_net

#k折交叉验证所对应的结果
best_net=printing_Kfold_scores(train_data,train_label)             
#  4.计算随机森林的准确率
print("训练集：",best_net.score(train_data,train_label))
print("测试集：",best_net.score(test_data,test_label))            
file = open(SavePath, "wb")
#将模型写入文件：
pickle.dump(best_net, file)#pickle.dumps()将对象obj对象序列化并返回一个byte对象
#最后关闭文件：
file.close()            

3、 模型测试

"""
import numpy as np
import gdal
import pickle 

#读取tif数据集
def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName+"文件无法打开")
    return dataset


#保存tif文件函数，首席判断数据类型，在对数据维度进行处理
def writeTiff(im_data,im_geotrans,im_proj,path):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    elif len(im_data.shape) == 2:
        im_data = np.array([im_data])
        im_bands, im_height, im_width = im_data.shape
    #创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(path, int(im_width), int(im_height), int(im_bands), datatype)
    if(dataset!= None):
        dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
        dataset.SetProjection(im_proj) #写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])
    del dataset
    
RFpath = r"E:\LZ\Data\model.pickle"
Landset_Path = r"E:\LZ\ROI\growth.tif"  #遥感图像
SavePath = r"E:\LZ\ROI\train_roi\train\save1.tif"
dataset = readTif(Landset_Path)#首先将栅格图像读取


Tif_width = dataset.RasterXSize #栅格矩阵的列数
Tif_height = dataset.RasterYSize #栅格矩阵的行数
Tif_bands = dataset.RasterCount #波段数
Tif_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息
Tif_proj = dataset.GetProjection()#获取投影信息
Landset_data = dataset.ReadAsArray(0,0,Tif_width,Tif_height)
print(Landset_data.shape)

##############################################调用保存好的模型
#以读二进制的方式打开预训练的模型文件
file = open(RFpath, "rb")
#把模型从文件中读取出来
svm_model = pickle.load(file)
#关闭文件
file.close()
################################################用读入的模型进行预测
#  在与测试前要调整一下数据的格式
data = np.zeros((Landset_data.shape[0],Landset_data.shape[1]*Landset_data.shape[2]))#波段数，行列
for i in range(Landset_data.shape[0]):
    data[i] = Landset_data[i].flatten() 

data = data.swapaxes(0,1)
print(data.shape)
data[np.isnan(data)] = 0
#for i in range(data.shape[0]):
#    for j in range(data.shape[1]):
#        if np.isnan(data[i][j]).any()=="True":
#            data[i][j]=0

#print(np.isnan(data).any())
#对调整好格式的数据进行预测
pred = svm_model.predict(data)
print(np.unique(pred))
#  同样地，我们对预测好的数据调整为我们图像的格式
pred = pred.reshape(Landset_data.shape[1],Landset_data.shape[2])*255
print(pred)
pred = pred.astype(np.uint8)

#  将结果写到tif图像里
writeTiff(pred,Tif_geotrans,Tif_proj,SavePath)