利用python进行遥感影像分类结果精度评估全流程

随夢而

一、引言

    最近利用python实现了遥感影像分类后结果的精度评估工作，切实感受到做算法开发时思路的重要性，有了好的思路，后期的编程工作如抽丝剥茧一般逐步实现，本次的精度评估方法主要思路如下：
    本文对提出的算法分类结果进行精度验证，但是缺少实测站点数据，因此利用人工视觉解译的结果作为标准分类结果。首先对GF-3遥感影像进行人工视觉解译，得到研究区海冰的分类图，基于格网的统计方法来验证算法的分类精度。具体为：利用ArcGIS在整个湖区建立格网，格网大小为1200×1200m2。分别计算各格网内算法分类结果和目视解译分类结果的面积，通过计算两者之间的相关系数、均方根误差等统计量，对分类结果进行精度评价，预期得到的结果类似于下图所示（原理相通，下图放了一个之前做过的案例）：



图1 算法分类结果与目视解译结果对比



图2 最终精度评估的散点图

二、输入数据

    由上可知本文要做的工作是对GF3 SAR数据提取的海冰结果进行精度验证，利用人工视觉解译的结果作为验证标准，代码集成思路如下：
    step1_创建渔网，基于算法得到的分类结果建立渔网；
    step2_首先基于外缘线建立缓冲区，然后利用渔网与外缘线缓冲区矢量数据进行叠加分析，提取一定缓冲区范围的渔网网格；
    step3_算法得到的分类矢量数据转栅格数据以及目视解译得到的矢量数据转栅格数据；
    step4_进行分区统计，利用step2得到的渔网矢量数据去统计step3得到的分类结果栅格数据；
    step5_画散点图进行展示。
    输入数据为：
    1、基于算法得到的分类结果矢量数据
    2、目视解译得到的矢量数据
    3、冰架外缘线矢量数据



图3 输入数据

三、python代码

    程序代码都已封装好且可以正常运行，但是有许多地方还需要优化和改进，今天写这篇博客的目的就是为了记录集成的思路和过程。代码运行前需要配置好运行环境，比较几个重要第三包（gdal、geopandas、rasterstats...）的安装需要特别注意，之前的博客中有记录，直接上代码，直接可以运行。
# -*- coding: utf-8 -*-
"""
@autho:YuanRuiQiang
@date:2022/9/6 16:11
@Location:BeiJing_PieSat
@breif:进行海冰提取的精度验证，指标要求精度优于90%
集成开始日期为20220906
step1_创建渔网，基于算法得到的分类结果建立渔网数据
step2_首先基于外缘线建立缓冲区，然后利用渔网与外缘线缓冲区矢量数据进行叠加分析，提取一定缓冲区范围的的渔网网格
step3_分类结果矢量转栅格以及目视解译矢量数据转栅格
step4_进行分区统计，利用step2得到的渔网去统计step3得到的分类结果栅格数据
step5_画散点图进行展示
"""
from osgeo import ogr, os, osr
import os
from math import ceil
import sys
from osgeo import gdal
import geopandas as gpd
import pandas as pd
import numpy as np
import rasterio
import rasterstats
import seaborn as sns
import warnings
import matplotlib as mpl
import matplotlib.pyplot as plt

def read_img(filename):
    # 思路主要是读取原始的GF3影像得到原始影像的范围信息和投影信息
    dataset = gdal.Open(filename)

    width = dataset.RasterXSize
    height = dataset.RasterYSize
    band = dataset.RasterCount
    im_data = dataset.ReadAsArray(0, 0, width, height)

    geotrans = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    return im_data, width, height, proj, geotrans

########################### step1创建渔网代码思路 ########################
# 1、输入渔网范围，格网宽度和高度
# 2、计算行列数
# 3、初始化左上角格网四角范围，创建输出文件
# 4、遍历每一列：
#         初始化第一行格网上下边界（放在行循环外，以便每次进入下一列定位到第一行）
#         遍历每一行：
#                 计算单个格网范围
#                 创建ring，写入格网四点，创建polygon，写入ring，构成一个格网
#                 多边形写入图层
#                 本行单个格网写入完成，行数加一并向下计算下一格网上下边界
#         本列格网全部写入，列数加一向右计算下一列格网左右边
# 创建渔网
def Fishgrid(outfile,xmin,xmax,ymin,ymax,gridwidth,gridheight):
    # 参数转换到浮点型
    xmin = float(xmin)
    xmax = float(xmax)
    ymin = float(ymin)
    ymax = float(ymax)
    gridwidth = float(gridwidth)
    gridheight = float(gridheight)
    # 计算行数和列数
    rows = ceil((ymax-ymin)/gridheight)
    cols = ceil((xmax-xmin)/gridwidth)
    # 初始化起始格网四角范围
    ringXleftOrigin = xmin
    ringXrightOrigin = xmin+gridwidth
    ringYtopOrigin = ymax
    ringYbottomOrigin = ymax-gridheight

    # 创建输出文件
    outdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outfile):
        outdriver.DeleteDataSource(outfile)
    outds = outdriver.CreateDataSource(outfile)
    # 辅助矢量数据信息读取！！！
    ds = ogr.Open(r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\worldboundary\GF3_KRN_NSC_029602_E71.5_S67.5_20220325_L2_HV_L2000636Polygon.shp')
    layer = ds.GetLayer(0)
    spatialRef = layer.GetSpatialRef()
    # print(spatialRef)
    extent = layer.GetExtent()
    print('extent:', extent)
    # print('upleft:', extent[0], extent[3])
    # print('bottomright:', extent[1], extent[2])

    outlayer = outds.CreateLayer(outfile, spatialRef, geom_type=ogr.wkbPolygon)  # 使用OGR 创建Shapefile,yrq添加了辅助数据的投影信息
    if outlayer == None:
        print("图层创建失败！")
        return
    # 不添加属性信息，获取图层属性
    outfielddefn = outlayer.GetLayerDefn()
    # 遍历列，每一列写入格网
    col = 0
    while col < cols:
        # 初始化，每一列写入完成都把上下范围初始化
        ringYtop = ringYtopOrigin
        ringYbottom = ringYbottomOrigin
        # 遍历行，对这一列每一行格子创建和写入
        row = 0
        while row < rows:
            # 创建左上角第一个格子
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYbottom)
            ring.CloseRings()
            # 写入几何多边形
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)
            # 创建要素，写入多边形
            outfeat = ogr.Feature(outfielddefn)
            outfeat.SetGeometry(poly)
            # 写入图层
            outlayer.CreateFeature(outfeat)
            outfeat = None
            # 下一多边形，更新上下范围
            row += 1
            ringYtop = ringYtop - gridheight
            ringYbottom = ringYbottom-gridheight
        # 一列写入完成后，下一列，更新左右范围
        col += 1
        ringXleftOrigin = ringXleftOrigin+gridwidth
        ringXrightOrigin = ringXrightOrigin+gridwidth
    # 写入后清除缓存
    outds = None

##################step2_渔网与外缘线缓冲区叠加分析#################
# 思路：先进行缓冲区分析，然后执行叠加分析
# inShp_path为输入shp的路径
# outShp_path为缓冲区分析后的输出shp的路径
# bufferDist表示缓冲区距离
def createBuffer(inShp_path, outShp_path, bufferDist):
    inputds = ogr.Open(inShp_path)
    inputlyr = inputds.GetLayer()
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outShp_path):
        driver.DeleteDataSource(outShp_path)
    ds = driver.CreateDataSource(outShp_path)
    layer = inputds.GetLayer(0)
    spatialRef = layer.GetSpatialRef()
    layer = ds.CreateLayer(outShp_path, spatialRef, geom_type=ogr.wkbPolygon)
    featureDefn = layer.GetLayerDefn()
    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)
        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        layer.CreateFeature(outFeature)
        outFeature = None

def clip_parking(clipfile_path,Data_path,OutPut_path):
    try:
        # 裁剪函数
        Maskfeture = gpd.read_file(clipfile_path)
        Indata = gpd.read_file(Data_path)
        Maskfeture = Maskfeture.to_crs(Indata.crs)
        Clipresult = gpd.clip(Indata, Maskfeture)
        # 保存文件SHP
        # save_SHP_path = OutPut_path + "\\" + "ClipResult.shp"
        Clipresult.to_file(OutPut_path, driver="ESRI Shapefile", encoding="utf-8")
    except:
        pass
#############################step3_矢量转栅格#######################
# 思路：将算法得出的结果和目视解译得到的结果矢量统一转为栅格数据进行下一步的分区统计
# 1.输入数据为算法检测完的矢量结果和目视解译的矢量结果
# 2.输出数据为上述两者的栅格图层
def ShapefileToRaster(vector_fn, proj, raster_fn):

    # Define pixel_size and NoData value of new raster
    pixel_size = 20
    NoData_value = 0

    # Open the data source and read in the extent
    source_ds = ogr.Open(vector_fn)
    source_layer = source_ds.GetLayer()
    # spatialRef = source_layer.GetSpatialRef()
    # print(spatialRef)
    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    # Create the destination data source
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Float32)
    target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
    target_ds.SetProjection(proj)

    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(NoData_value)
    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])

##########################step4_分区统计_面积制表######################
# 1.利用rasterstats包进行分区统计，rasterstats包的使命就是基于矢量数据进行栅格数据的统计分析
# 2.将计算后的面积作为返回值返回

def Zonal_Statistics(districts,raster1,raster2):
    try:
        ice_data1 = raster1.read(1)
        # 设置坐标变换信息
        affine1 = raster1.transform
        # 准备开始进行空间分区计算
        # 第一个参数是矢量分区，第二个是栅格，第三个是坐标变换信息，第四个是统计均值
        sta_sum1 = rasterstats.zonal_stats(districts, ice_data1, affine=affine1, stats=['sum'], geojson_out=True)
        ice_data2 = raster2.read(1)
        # 设置坐标变换信息
        affine2 = raster2.transform
        # 准备开始进行空间分区计算
        # 第一个参数是矢量分区，第二个是栅格，第三个是坐标变换信息，第四个是统计均值
        sta_sum2 = rasterstats.zonal_stats(districts, ice_data2, affine=affine2, stats=['sum'], geojson_out=True)
        test1 = [sta_sum1['properties']['sum'] for i in range(0, len(sta_sum1))]
        test2 = [sta_sum2['properties']['sum'] for i in range(0, len(sta_sum2))]
        return test1, test2
    except:
        pass
############################step5_画散点图进行展示######################
# 1.将上个程序计算后的面积的返回值作为输入值
# 2.绘制散点图
# 3.绘制趋势线，添加R方
def Curve_Fitting(x,y):
    plt.scatter(x, y, c='blue', alpha=0.5, marker='*')
    plt.xlabel('GF3_thredhold', c='black')
    plt.ylabel('visual', c='black')
    plt.xlim((0, 800))
    plt.ylim((0, 800))
    plt.tick_params(axis='x', colors='black')
    plt.tick_params(axis='y', colors='black')
    plt.show()

if __name__ == '__main__':
    print("start")
    os.chdir(r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2')
    shp = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\fishgrid.shp'
    outline = 'GF3_KRN_NSC_029602_polyline.shp'
    outShp_path = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\buffer\butter_polygon.shp'
    outFile_path = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\clip\ClipResult.shp'
    bufferDist = 1500
    print("start process fishnet")
    Fishgrid(shp, 3744316, 4207996, 735825, 1133525, 1200, 1200)  # 执行生成渔网
    createBuffer(outline, outShp_path, bufferDist)  # 执行缓冲区分析
    print("start process clip")
    clip_parking(outShp_path, shp, outFile_path)  # 进行空间叠加分析
    print("start process shpToRaster")
    raster = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\1_GF3_extractResult.tif'  # 引入原始影像的投影
    im_data, width, height, proj, geotrans = read_img(raster)
    vector_fn1 = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\ICE_result\GF3_threshold.shp'
    vector_fn2 = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\ICE_result\Visual_interpretation.shp'
    raster_fn1 = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\vectorToRaster\test_1.tif'
    raster_fn2 = r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\vectorToRaster\test_2.tif'
    ShapefileToRaster(vector_fn1, proj, raster_fn1)
    ShapefileToRaster(vector_fn2, proj, raster_fn2)
    # 使用geopandas读取矢量数据
    districts = gpd.read_file(r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\clip\ClipResult.shp')
    # 使用rasterio读取栅格数据，栅格数据和矢量数据的坐标投影需要一致
    print("start Zonal_Statistics")
    raster1 = rasterio.open(r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\vectorToRaster\test_1.tif')
    raster2 = rasterio.open(r'G:\Work\projects\2022-GF2-ICE\ice_integtration3\Validation\part2\vectorToRaster\test_2.tif')
    test1, test2 = Zonal_Statistics(districts, raster1, raster2)
    Curve_Fitting(test1, test2)

    print("All done")四、输出结果




图4 最终精度评估散点图

庄洪涛

好帖必须得顶起

彩色风

不错 支持下

空虚

大人，此事必有蹊跷！