利用Python的Scikit-Learn對遙感影像進行隨機森林分類

2021-03-02 開放GIS實驗室與防災減災服務

隨機森林是一個包含多個決策樹的分類器,因其運算速度快、分類精度高、算法穩定等特點,被廣泛應用到遙感圖像的分類研究中。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]

                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]

                X=leftX+j*rX+.00001

                if point1.X==point2.X:

                    intersectX=point1.X

                    if intersectX>X:count+=1

                else:

                   k=(point2.Y-point1.Y)/(point2.X-point1.X)

                    if k==0:

                        if X<point1.X or X<point2.X:

                            count+=1

                    else:

                       intersectX=(Y-point1.Y)/k+point1.X

                        ifintersectX>X:count+=1

 

            if count%2==0:

                if BandsCount>1:

                    for bc in range(BandsCount):

                       arr[bc][i][j]=(nodatavalue)

                else:

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

   rasterspatial = 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 False

   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=False

                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

五、隨機森林分類

from sklearn.ensemble 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]#左上角點

   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 = tbrf.createClassifier(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 = math.ceil(height / bl)

       if blockSize < 1:blockSize = 1

       if blockSize > 1:blockSize+=5

   if  blockSize != 1:

       blockHeight = 0

       modHeight = 0

       modHeight = height % blockSize

 

       if modHeight == 0:

            blockHeight = int(height / blockSize)

           

       else:

            blockHeight = int(height / blockSize)

       print(f"分塊大小{width}*{blockHeight}")

       for bs in range(blockSize):

            print(f"計算塊{bs+1}/{blockSize}")

           checkMemory(500)

            arr = rds.ReadAsArray(0,bs * blockHeight,width,blockHeight)

           

            for i in range(blockHeight):

                print(f"分塊:{bs+1}/{blockSize}添加分類數據{round((i+1)*100/blockHeight,4)}%")

                for k in range(width):

                    tem = list()

                    for bc inrange(BandsCount):

                        tem.append(int(arr[bc][i][k]))

                    fixX.append(tem)

            print(f"分塊:{bs+1}/{blockSize}計算分類結果……")

            checkMemory(800)

            z = clf.predict(fixX)        

            Z.extend(z)

            fixX = list()

            arr = None

       print(f"計算餘數:{width}*{modHeight}")

       checkMemory(500)

       arr = rds.ReadAsArray(0,blockSize *blockHeight,width,modHeight)

       if modHeight > 0:

            for i in range(modHeight):

                print(f"餘塊:添加分類數據{round((i+1)*100/modHeight,4)}%")

                for k in range(width):

                    tem = list()

                   for bc in range(BandsCount):

                        tem.append(int(arr[bc][i][k]))

                    fixX.append(tem)

            print("餘塊:計算分類結果……")

            checkMemory(500)

            z = clf.predict(fixX)        

            Z.extend(z)

       Z = numpy.array(Z)

       #Z=Z.reshape(1,width*height)

       Z = Z.reshape(height,width)

       fixX = None

       arr = None

   else:

       checkMemory(1000)

       arr = rds.ReadAsArray(0,0,width,height)

       for i in range(height):

            print(f"添加訓練樣本{round((i+1)*100/height,4)}%")

            for k in range(width):

                tem = list()

                for bc in range(BandsCount):

                    tem.append(int(arr[bc][i][k]))

                fixX.append(tem)

       arr = None

       print("計算分類結果……")

       Z = clf.predict(fixX)

       Z = numpy.array(Z)

       Z = Z.reshape(height,width)

   driver = gdal.GetDriverByName("GTiff")

   filepath,filename = os.path.split(outfile)

   short,ext = os.path.splitext(filename)

   print("創建輸出文件")

   out = driver.Create(outfile,width,height,1,rds.GetRasterBand(1).DataType)

   out.SetGeoTransform(transform)

   out.SetProjection(rds.GetProjectionRef())

   print("寫入數據……")

   out.GetRasterBand(1).WriteArray(Z)

   out.FlushCache()

   out = None

print("計算完成")

六、分類圖像和結果

點擊閱讀原文查看數據獲取方法    

相關焦點

  • Scikit-learn估計器分類
    Scikit-learn實現了一系列數據挖掘算法,提供通用編程接口、標準化的測試和調參公局,便於用戶嘗試不同算法對其進行充分測試和查找最優參數值。本次講解數據挖掘通用框架的搭建方法。有了這樣一個框架,後續章節就可以把講解重點放到數據挖掘應用和技術上面。為幫助用戶實現大量分類算法,scikit-learn把相關功能封裝成估計器。估計器用於分類任務,它主要包括兩個函數。
  • WePay機器學習反欺詐實踐:Python+scikit-learn+隨機森林
    WePay採用了流行的Python、scikit-learn開源學習機器學習工具以及隨機森林算法。以下是文章內容:什麼是shellselling?雖然欺詐幾乎涉及各種領域,但相對於傳統的買方或賣方僅僅擔心對方是否是騙子,支付平臺需要擔心的是交易雙方。如果其中任何一方存在信用詐騙,真正的持卡人發現和撤銷費用,平臺自身就要進行帳單償還。
  • 使用Scikit-learn 理解隨機森林
    scikit-learn,作者 ando。我的一些代碼包正在做相關工作,然而,大多數隨機森林算法包(包括 scikit-learn)並沒有給出預測過程的樹路徑。因此 sklearn 的應用需要一個補丁來展現這些路徑。
  • 第二章 使用scikit-learn進行分類預測
    01scikit-learn estimatorsestimators,允許算法的標準化實現和測試一個通用的、輕量級的接口,供分類器遵循。通過使用這個接口,我們可以將這些工具應用於任意分類器,而不需要擔心算法如何工作。
  • Python中的人工智慧入門:在scikit-learn中的建模
    問題是根據萼片和花瓣寬度和長度的測量結果對三種鳶尾花中的一種進行分類。此外,我們將使用seaborn和scikit-plot庫進行可視化,所以我們也將安裝它們。!pip install seaborn scikit-plot使用如此知名的數據集的方便之處在於,我們可以很容易地從很多包中加載數據集,例如,像這樣。
  • 一遍就懂的Scikit-Learn機器學習分類過程
    為了對鳶尾花進行分類,我們將執行以下步驟:加載鳶尾花數據集將數據拆分為訓練集和測試集在訓練集上訓練模型計算模型精度對外部輸入數據執行預測為了執行所有這些步驟,我們將使用scikit-learn python庫:https://scikit-learn.orgscikit-learn
  • 利用Spark 和 scikit-learn 將你的模型訓練加快 100 倍
    當我們使用 Spark 進行數據處理時,我們首選的機器學習框架是 scikit-learn。隨著計算機變得越來越便宜,機器學習解決方案的上市時間變得越來越關鍵,我們探索了加快模型訓練的各種方法。其中一個解決方案是將 Spark 和 scikit-learn 中的元素組合到我們自己的混合解決方案中。
  • python機器學習之使用scikit-learn庫
    引言數據分析由一連串的步驟組成,對於其中預測模型的創建和驗證這一步,我們使用scikit-learn這個功能強大的庫來完成。scikit-learning庫python庫scikit-learn整合了多種機器學習算法。
  • 使用 scikit-learn 玩轉機器學習——集成學習
    Random Forests(隨機森林)我們都知道森林是由樹構成的(手動滑稽,QAQ),所以隨機森林也不例外,隨機森林裡面的樹叫做決策樹。上次我們剛聊過決策樹,相信小夥伴們還有些印象,決策樹是由一系列節點構成的,每劃分一個節點都要在所有的特徵維度的每個特徵可能取到的值上進行搜索,以取得信息熵的最小和,或最大的信息增益。
  • 隨機森林入門攻略
    Amit, Gemen和Ho Tim Kam各自獨立地介紹了特徵隨即選擇的思想,並且運用了Breiman的「套袋」思想構建了控制方差的決策樹集合。在此之後,Deitterich在模型中引入了隨即節點優化的思想,對隨機森裡進行了進一步完善。什麼是隨機森林?隨機森林是一種多功能的機器學習算法,能夠執行回歸和分類的任務。
  • Scikit-learn使用總結
    在機器學習和數據挖掘的應用中,scikit-learn是一個功能強大的python包。在數據量不是過大的情況下,可以解決大部分問題。學習使用scikit-learn的過程中,我自己也在補充著機器學習和數據挖掘的知識。這裡根據自己學習sklearn的經驗,我做一個總結的筆記。另外,我也想把這篇筆記一直更新下去。
  • Python粉都應該知道的開源機器學習框架:Scikit-learn入門指南
    和其他眾多的開源項目一樣,Scikit-learn目前主要由社區成員自發進行維護。可能是由於維護成本的限制,Scikit-learn相比其他項目要顯得更為保守。這主要體現在兩個方面:一是Scikit-learn從來不做除機器學習領域之外的其他擴展,二是Scikit-learn從來不採用未經廣泛驗證的算法。
  • 極簡Scikit-Learn入門
    Scikit-Learn高清全景圖傳送:http://scikit-learn.org/stable/tutorial/machine_learning_map/index.html在機器學習和數據挖掘的應用中
  • 用Python (scikit-learn) 做PCA分析
    Iris數據集是scikit-learn附帶的數據集之一,不需要從外部網站下載任何文件。下面的代碼將加載iris數據集。如果你想看到不縮放數據可能帶來的負面影響,scikit-learn有一節是講關於不標準化數據的影響的(https://scikit-learn.org/stable/auto_examples/preprocessing/plot_scaling_importance.html#sphx-glr-auto-examples-preprocessing-plot-scaling-importance-py
  • python實現隨機森林
    定義:隨機森林指的是利用多棵決策樹對樣本進行訓練並預測的一種分類器。可回歸可分類。
  • Python + Scikit-learn 完美入門機器學習指南 ​
    如果你是一名 Python 程式設計師,又正好想學習一下人工智慧技術,scikit-learn 可能是你最好的選擇之一。實驗樓上線了一門新課 —《scikit-learn 機器學習入門實戰》,可以帶大家快速掌握 scikit-learn 框架使用方法、機器學習基礎概念、常用算法等知識點,並通過模型的搭建和選擇,進一步提升你的實踐能力。
  • scikit-learn機器學習簡介
    監督學習,數據包含有我們要預測的額外屬性(單擊此處 轉到scikit-learn監督學習頁面)。其中包括:分類:樣本屬於兩個或多個類別,我們想從已經標記的數據中學習如何預測未標記數據的類別。分類問題的一個例子是手寫數字識別,手寫數字識別的目的是將每個輸入向量分配給有限的離散類別中的一個。
  • 開源機器學習框架:Scikit-learn API簡介
    分類:是指識別給定對象的所屬類別,屬於監督學習的範疇,最常見的應用場景包括垃圾郵件檢測和圖像識別等。目前Scikit-learn已經實現的算法包括:支持向量機(SVM),最近鄰,邏輯回歸,隨機森林,決策樹以及多層感知器(MLP)神經網絡等等。
  • 教程 | 如何在Python中用scikit-learn生成測試數據集
    scikit-learn Python 庫提供一套函數,用於從可配置測試問題中生成樣本來進行回歸和分類。在本教程中,你將學習測試問題及如何在 Python 中使用 scikit-learn 進行測試。scikit-learn 是一個用於機器學習的 Python 庫,它提供了生成一組測試問題的函數。在本教程中,我們將看一些為分類和回歸算法生成測試問題的例子。分類測試問題分類是將標籤分配給數據的問題。
  • 利用 Scikit Learn的Python數據預處理實戰指南
    簡而言之,預處理是指在你將數據「餵給」算法之前進行的一系列轉換操作。在Python中,scikit-learn庫在sklearn.preprocessing下有預裝的功能。有更多的選擇來進行預處理,這將是我們要探索的。讀完本文,你將具備數據預處理的基本技能並對其有更深入的理解。