python随机森林_利⽤Python的ScikitLearn对遥感影像进⾏随
机森林分类
随机森林是⼀个包含多个决策树的分类器，因其运算速度快、分类精度⾼、算法稳定等特点，被⼴泛应⽤到遥感图像的分类研究中。Scikit-Learn作为Python 编程语⾔的免费软件机器学习库，提供了对随机森林算法的⽀持，但没有提供针对遥感影像分类的相关函数。因此，本篇⽂章将为读者介绍利⽤Python及其扩展包Scikit-Learn对遥感影像进⾏随机森林分类的完整过程，包括：ShapeFile格式样本数据的读取、栅格数据读取和裁剪、利⽤Scikit-Learn的RandomForestClassifier模块进⾏样本训练和遥感影像分类。
⼀、Scikit-Learn的安装
直接执⾏命令pip install scikit-learn，所有依赖库都会⾃动安装。安装完成后，添加代码
ble import RandomForestClassifier即可使⽤。
⼆、样本绘制
在ArcGIS中绘制训练样本，格式为shpfile，可以是点类型或线类型，建⽴Value字段，⽤于存储分类编号，例如1-湿地，2-湖泊，3-⽔稻。
三、栅格数据裁剪
通过shp样本获取对应的栅格值，需要使⽤多边形裁剪栅格数据，我们使⽤射线算法。全部代码如下(TrainByRandomForest.py)：
ble 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
if point.X>maxX:maxX=point.X
宁波一对一
if point.Y
if point.Y>maxY:maxY=point.Y
leftX=minX
upY=maxY
shy什么意思
rightX=maxX
martinibottomY=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]
point2=polygonPoints[k+1]
if (point1.Y>=Y and point2.Y<=Y) or (point1.Y<=Y and point2.Y>=Y): pointsindex.append(k)
for j in range(width):
count=0
for m in (pointsindex):
point1=polygonPoints[m]
point2=polygonPoints[m+1]
you tubeX=leftX+j*rX+.00001
if point1.X==point2.X:
interctX=point1.X
if interctX>X:count+=1
el:
k=(point2.Y-point1.Y)/(point2.X-point1.X)
if k==0:
if Xor X
count+=1
el:
interctX=(Y-point1.Y)/k+point1.X
ifinterctX>X:count+=1
2017四六级成绩查询
if count%2==0:
if BandsCount>1:
for bc in range(BandsCount):
arr[bc][i][j]=(nodatavalue)
el:
arr[i][j]=-1
#为了测试结果的正确性，可以先将其写到硬盘
#WriteRaster("test.tif",arr,inraster,width,height,BandsCount,leftX,upY) return arr,width,height,BandsCount,leftX,upY
四、创建分类器并训练样本
def createClassifier(inraster,inshp,field:str="Id",treenum:int=100):
heal the worldrasterspatial = gdal.Open(inraster)
spatial2=osr.SpatialReference()
spatial2.ImportFromWkt(rasterspatial.GetProjectionRef())
shpspatial=ogr.Open(inshp)
layer=shpspatial.GetLayer(0)
spatial1=layer.GetSpatialRef()
ct=osr.CreateCoordinateTransformation(spatial1,spatial2)
oFeature = layer.GetNextFeature()
# 下⾯开始遍历图层中的要素
geom=oFeature.GetGeometryRef()
if geom.GetGeometryType()==ogr.wkbPoint:
return createClassifierByPoint(inraster,inshp)
k=geom.GetGeometryType()
if geom.GetGeometryType()!=ogr.wkbPolygon:
print("样本必须为单部件多边形")
孤儿院英文
return Fal
trainX = list()
trainY = list()
print("读取样本")
while oFeature is not None:
花式调酒培训
geom=oFeature.GetGeometryRef()
wkt=geom.ExportToWkt()
points=WKTToPoints(wkt)
polygonPoints=[]
value=oFeature.GetField(field)
for point in points:
pC=ct.TransformPoint(point.X,point.Y,0)
polygonPoints.append(Point(pC[0],pC[1]))
arr,width,height,BandsCount,leftX,upY=GetSubRaster(inraster,polygonPoints) for i in range(height):
for k in range(width):
nodata=True
tem = list()
for bc in range(BandsCount):
v=int(arr[bc][i][k])
tem.append(v)
if v>0:nodata=Fal
if nodata:
continue
trainX.append(tem)
trainY.append(int(value))
oFeature = layer.GetNextFeature()
ct=None
spatial1=None
spatial2=None
print("训练样本")
我最喜欢的一首诗clf = RandomForestClassifier(n_estimators=treenum)
clf.fit(trainX, trainY)#训练样本
print("训练完成")
return clf
五、随机森林分类
ble import RandomForestClassifier
import osgeo
from osgeo import gdal
from osgeo import osr
from osgeo import ogr
import numpy
import os
import sys
import TrainByRandomForest as tbrf
def RandomForestClassification(ClassifyRaster,TrainRaster,TrainShp,outfile,blockSize=0,treenum=100,max_depth=10): rds = gdal.Open(ClassifyRaster)
#print((rds.GetRasterBand(1).DataType))
transform = (rds.GetGeoTransform())
lX = transform[0]#左上⾓点den
lY = transform[3]
rX = transform[1]#分辨率
rY = transform[5]
width = rds.RasterXSize
height = rds.RasterYSize
bX = lX + rX * width#右下⾓点
bY = lY + rY * height
BandsCount = rds.RasterCount
clf = ateClassifier(TrainRaster,TrainShp)
Z = list()
fixX = list()
if blockSize == 0:
p,a = memory_usage()
pv = 0.6 / 10000
checkMemory(2000)#内存⼩于2GB，不在计算
bl = (a - 2000) / pv / height / BandsCount
blockSize = il(height / bl)
