利用Python的Scikit-Learn库对遥感影像进行随机森林分类

利用Python的Scikit-Learn库对遥感影像进行随机森林(RandomForest)分类【附视频教程】

**

视频教程地址：https://www.bilibili.com/video/BV17U4y187J4/

**

随机森林是一个包含多个决策树的分类器，因其运算速度快、分类精度高、算法稳定等特点，被广泛应用到遥感图像的分类研究中。Scikit-Learn作为Python 编程语言的免费软件机器学习库，提供了对随机森林算法的支持，但没有提供针对遥感影像分类的相关函数。因此，本篇文章将为读者介绍利用Python及其扩展包Scikit-Learn对遥感影像进行随机森林分类的完整过程，包括：ShapeFile格式样本数据的读取、栅格数据读取和裁剪、利用Scikit-Learn的RandomForestClassifier模块进行样本训练和遥感影像分类。

一、Scikit-Learn的安装

直接执行命令pip install scikit-learn，所有依赖库都会自动安装。安装完成后，添加代码from sklearn.ensemble import RandomForestClassifier即可使用

二、样本绘制

在ArcGIS中绘制训练样本，格式为shpfile，可以是点类型或面类型，建立Value字段，用于存储分类编号，例如1-湿地，2-湖泊，3-水稻。

三、栅格数据裁剪

通过shp样本获取对应的栅格值，需要使用多边形裁剪栅格数据，我们使用射线算法。全部代码如下（TrainByRandomForest.py）：

代码如下（示例）：

from sklearn.ensemble import RandomForestClassifier
import osgeo
from osgeo import gdal
from osgeo import osr
from osgeo import ogr
def GetSubRaster(inraster,polygonPoints:list):
    polygonPoints.append(polygonPoints[0])#面多边形坐标封闭
	print("当前多边形节点数量："+str(len(polygonPoints)))
	#计算最小边界矩形
    minX=10000000000000
    maxX=-minX
    minY=100000000000000000
    maxY=-minY

    for point in polygonPoints:
        if point.X<minX:minX=point.X
        if point.X>maxX:maxX=point.X
        if point.Y<minY:minY=point.Y
        if point.Y>maxY:maxY=point.Y
    leftX=minX
    upY=maxY
    rightX=maxX
    bottomY=minY

    rds = gdal.Open(inraster) 
    transform = (rds.GetGeoTransform())
    lX = transform[0]#左上角点
    lY = transform[3]
    rX = transform[1]#分辨率
    rY = transform[5]

    wpos=int((leftX-lX)/rX)
    hpos=int((upY-lY)/rY)

    width=int((rightX-leftX)/rX)
    height=int((bottomY-upY)/rY)
    BandsCount = rds.RasterCount
    arr = rds.ReadAsArray(wpos,hpos,width,height)
    fixX=list()
    nodatavalue=rds.GetRasterBand(1).GetNoDataValue()
    for i in range(height):
        if height>200:
            print(f"多边形裁剪进度：{round(((i+1)/height)*100,4)}%")
        Y=upY+i*rY+.00001
	    #射线算法只需要比对多边形的一条水平线上的边
        pointsindex=list()
        for k in range(len(polygonPoints)-1):
                point1=polygonPoints[k]</