随机森林遥感影像识别+Rater转MultiShp

目录

随机森林代码参考，在此表示感谢

代码运行遇到的问题：

一.有些遥感影像虽然是tif格式的，但是直接运行会报错，把遥感影像在ArcGIS中重新导出数据，适用于国产数据

二.osgeo的gdal包的安装问题：

1、不要在anaconda中直接安装GDAL，先下一个XXX.whl在安装，参考上面的链接

2、Cannot find proj.db的错误

3、用python打包成exe的时候，python 3.9 + gdal 3.3.2版本会报错，回退到python3.7 + gdal 3.3.1就好了（提醒打包失败，一般都是环境的问题）

Rater转MultiShp

小图斑去除

栅格矢量化

按分类结果提取每一类的矢量文件

矢量平滑

异步问题！！！

感谢

---

随机森林代码参考，在此表示感谢

利用Python的Scikit-Learn库对遥感影像进行随机森林分类_u014444411的博客-CSDN博客利用Python的Scikit-Learn库对遥感影像进行随机森林(RandomForest)分类随机森林是一个包含多个决策树的分类器，因其运算速度快、分类精度高、算法稳定等特点，被广泛应用到遥感图像的分类研究中。Scikit-Learn作为Python 编程语言的免费软件机器学习库，提供了对随机森林算法的支持，但没有提供针对遥感影像分类的相关函数。因此，本篇文章将为读者介绍利用Python及其扩展包Scikit-Learn对遥感影像进行随机森林分类的完整过程，包括：ShapeFile格式样本数据的读取、https://blog.csdn.net/u014444411/article/details/108357109?spm=1001.2014.3001.5501

代码运行遇到的问题：

一.有些遥感影像虽然是tif格式的，但是直接运行会报错，把遥感影像在ArcGIS中重新导出数据，适用于国产数据

二.osgeo的gdal包的安装问题：

anaconda3安装GDAL - 知乎一、GDAL (Geospatial Data Abstraction Library) GDAL库最初是用C和C ++编写的，但它与其他几种语言，包括Python，做了绑定，所以尽管这些代码没有用Python进行重写，但它为Python中使用GDAL/OGR库提供了接口。因…https://zhuanlan.zhihu.com/p/429789848

1、不要在anaconda中直接安装GDAL，先下一个XXX.whl在安装，参考上面的链接

2、Cannot find proj.db的错误

ERROR 1: PROJ: proj_identify: Cannot find proj.db
ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db

 解决方法一：把GDAL中的proj文件添加到系统变量中

解决方法二：

如果打包成exe换个电脑，就又会缺少环境变量，又会报这个错误，

解决办法，把proj文件复制到exe同一文件下，在代码中添加如下代码：

path1 = sys.executable
path2 = os.path.dirname(path1)
path = os.path.join(path2, 'proj')
os.environ['PROJ_LIB'] = path
print(os.environ['PROJ_LIB'])

 

3、用python打包成exe的时候，python 3.9 + gdal 3.3.2版本会报错，回退到python3.7 + gdal 3.3.1就好了（提醒打包失败，一般都是环境的问题）

Rater转MultiShp

小图斑去除

其实做一个闭运算更好，但是GDAL没提供具体的算法，只能用cv2的算法，请自行百度把

from osgeo import gdal
 

def callback(v1, v2, v3):
    print(v1)


# Remove small spots
def RemoveSpots(inRaster, threshold=20, connectedness=4):
    print(inRaster)
    inputFile = inRaster
    options = ()
    mask = 'none'
    src_ds = gdal.Open(inputFile, gdal.GA_Update)
    srcBand = src_ds.GetRasterBand(1)
    dstBand = srcBand
    if mask is 'default':
        maskBand = srcBand.GetMaskBand()
    elif mask is 'none':
        maskBand = None
    else:
        mask_ds = gdal.Open(mask)
        maskBand = mask_ds.GetRasterBand(1)

    # 参数说明 输入数据波段 、设置掩码波段（只对掩码区域进行处理）、输出数据波段、去除板块的最大像元个数、图斑连通方式、占位(options)、进度条
    result = gdal.SieveFilter(srcBand, maskBand, dstBand, threshold, connectedness, options, callback)
    return result


if __name__ == '__main__':
    gdal.AllRegister()
    inRaster = "C:/Users/z6q6k6/Desktop/aa/result_1.tif"
    RemoveSpots(inRaster)

栅格矢量化

from osgeo import gdal
from osgeo import ogr
from pathlib import Path
from osgeo import osr
import os
import sys


def callback(v1, v2, v3):
    print(v1)


def RasterToShp(inputFile, dst_filename):
    driver = ogr.GetDriverByName('ESRI Shapefile')

    if os.path.exists(dst_filename):
        driver.DeleteDataSource(dst_filename)

    ds = gdal.Open(inputFile, gdal.GA_ReadOnly)
    srcband = ds.GetRasterBand(1)
    maskband = srcband.GetMaskBand()
    prj = osr.SpatialReference()
    prj.ImportFromWkt(ds.GetProjection())
    drv = ogr.GetDriverByName('ESRI Shapefile')
    dst_ds = drv.CreateDataSource(dst_filename)
    dst_layername = 'out'
    dst_layer = dst_ds.CreateLayer(dst_layername, prj, ogr.wkbMultiPolygon)
    dst_fieldname = 'Value'
    fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
    dst_layer.CreateField(fd)
    dst_field = 0

    options = []
    # 参数  输入栅格图像波段\掩码图像波段、矢量化后的矢量图层、需要将DN值写入矢量字段的索引、算法选项、进度条回调函数、进度条参数
    print("栅格转shp完成")
    return gdal.Polygonize(srcband, maskband, dst_layer, dst_field, options, callback)


if __name__ == '__main__':
    gdal.AllRegister()
    inputFile = "C:/Users/z6q6k6/Desktop/aa/result_1.tif"
    dst_filename = "C:/Users/z6q6k6/Desktop/aa/allShp.shp"
    RasterToShp(inputFile, dst_filename)

按分类结果提取每一类的矢量文件

把包含的所有类别的shp复制多份，然后在每一份中，删除不符合要求的类别，就得到了每一类的矢量文件

这个方法目前看起来比较low，但是还是学到了很多东西，值得分享！

在这再提供一种思路：

getValue函数中，feature.GetFieldAsString('Value')获取类别

然后参考CopyShp函数，分别导入到新创建的矢量中

from osgeo import gdal
from osgeo import ogr
from pathlib import Path
from osgeo import osr
import os
import sys
import shutil

def CopyShp(inputfile, pathfile, str):
    os.chdir(pathfile)
    # 设置driver
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 打开输入的矢量
    inDs = driver.Open(inputfile, 1)
    if inDs is None:
        print("Could not open", inputfile)
        sys.exit(1)
    # 打开输入矢量的图册
    inLayer = inDs.GetLayer()

    # 检查所需创建的矢量是否已存在
    if os.path.exists(str):
        driver.DeleteDataSource(str)
    # 创建矢量
    outDs = driver.CreateDataSource(str)
    if outDs is None:
        print("Could not create file")
        sys.exit(1)
    # 创建图册
    prj = osr.SpatialReference()
    a = inLayer.GetSpatialRef()
    s = a.ExportToWkt()
    prj.ImportFromWkt(s)
    outLayer = outDs.CreateLayer(str, srs=prj, geom_type=ogr.wkbPolygon)

    # for inFeature in inLayer:
    #     outLayer.CreateFeature(inFeature)
    # 将输入矢量的属性应用到新矢量
    fieldDefn = inLayer.GetFeature(0).GetFieldDefnRef('Value')
    outLayer.CreateField(fieldDefn)
    featureDefn = outLayer.GetLayerDefn()

    inFeature = inLayer.GetNextFeature()
    while inFeature:
        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(inFeature.GetGeometryRef())
        for i in range(inFeature.GetFieldCount()):
            value = inFeature.GetField(i)
            outFeature.SetField(i, value)
        outLayer.CreateFeature(outFeature)
        inFeature.Destroy()
        outFeature.Destroy()
        inFeature = inLayer.GetNextFeature()

    inDs.Destroy()
    outDs.Destroy()
    return os.path.join(pathfile, str)


def repack(shpFileName, strValue):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    pFeatureDataset = driver.Open(shpFileName, 1)

    if pFeatureDataset is None:
        print("Could not open", pFeatureDataset)
        sys.exit(1)
    pFeaturelayer = pFeatureDataset.GetLayer()

    # strValue = 0
    strFilter = "Value != '" + str(strValue) + "'"
    pFeaturelayer.SetAttributeFilter(strFilter)

    pFeatureDef = pFeaturelayer.GetLayerDefn()
    pLayerName = pFeaturelayer.GetName()
    pFieldName = "Value"
    # pFieldIndex = pFeatureDef.GetFieldIndex(pFieldName)
    pFeatureDef.GetFieldIndex(pFieldName)
    for pFeature in pFeaturelayer:
        pFeatureFID = pFeature.GetFID()
        pFeaturelayer.DeleteFeature(int(pFeatureFID))
    strSQL = "REPACK " + str(pFeaturelayer.GetName())
    pFeatureDataset.ExecuteSQL(strSQL, None, "")
    pFeatureLayer = None
    pFeatureDataset = None
    return 0


def getValue(inShp):
    gdal.SetConfigOption("SHAPE_ENCODING", "")
    ogr.UseExceptions()
    in_ds = ogr.Open(inShp)
    in_lyr = in_ds.GetLayer()
    list = []
    for feature in in_lyr:
        # print(feature.SetAttributeFilter("Value = '1'"))
        if not feature.GetFieldAsString('Value') in list:
            list.append(feature.GetFieldAsString('Value'))
    del in_ds
    return list


def SimpleToMulti(inputfile, pathfile):
    if os.path.exists(pathfile):
        # os.remove(outfile)    删除的时非空文件夹，os.remove会出现拒绝访问
        shutil.rmtree(pathfile)

    os.makedirs(pathfile)

    list = getValue(inputfile)
    sum = len(list)
    num = 0
    for i in list:
        pathAndName = CopyShp(inputfile, pathfile, i + '.shp')
        repack(pathAndName, i)
        num = num + 1
        print(f"处理完成{round(num / sum, 2) * 100}%")


if __name__ == '__main__':
    gdal.AllRegister()
    inputfile = "C:/Users/z6q6k6/Desktop/aa/allShp.shp"
    pathfile = "C:/Users/z6q6k6/Desktop/aa//eachShp"
    SimpleToMulti(inputfile, pathfile)

矢量平滑

用两个缓冲区（正缓冲区和负缓冲区）来平滑矢量的线

但是有一个问题，就是如果geometry.GetArea()太大，速度会很慢

from osgeo import gdal
from osgeo import ogr
from pathlib import Path
import tqdm
from osgeo import osr
import os
import sys


def callback(v1, v2, v3):
    print(v1)


def smoothing(inShp, fname, bdistance=0.001):
    """
    :param inShp: 输入的矢量路径
    :param fname: 输出的矢量路径
    :param bdistance: 缓冲区距离
    :return:
    """
    ogr.UseExceptions()
    in_ds = ogr.Open(inShp)
    in_lyr = in_ds.GetLayer()
    # 创建输出Buffer文件
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if Path(fname).exists():
        driver.DeleteDataSource(fname)
    # 新建DataSource，Layer
    out_ds = driver.CreateDataSource(fname)
    out_lyr = out_ds.CreateLayer(fname, in_lyr.GetSpatialRef(), ogr.wkbPolygon)
    def_feature = out_lyr.GetLayerDefn()
    # 遍历原始的Shapefile文件给每个Geometry做Buffer操作
    for feature in tqdm.tqdm(in_lyr):
        geometry = feature.GetGeometryRef()
        # print(geometry.GetArea())
        buffer = geometry.Buffer(bdistance).Buffer(-bdistance)
        out_feature = ogr.Feature(def_feature)
        out_feature.SetGeometry(buffer)
        out_lyr.CreateFeature(out_feature)
        out_feature = None
    out_ds.FlushCache()
    del in_ds, out_ds


if __name__ == '__main__':
    gdal.AllRegister()
    inputshp = "C:/Users/z6q6k6/Desktop/aa/eachShp/2.shp"
    output_path = "C:/Users/z6q6k6/Desktop/aa/eachShp/2_1.shp"
    smoothing(inputshp, output_path)
    print("end!")

异步问题！！！

在读取文件时候多是异步操作，如果上一步的文件还未完全创建（可能只有几k），就立即进行下一步的操作，一定会报错

感谢

本文参考了很多博客，不能一一添加引用，在此一并表示感谢