基于arcgis的遥感影像标签制作（目标检测）

文章目录

• ​​1. 在arcgis中新建线矢量​​
• ​​2. 检测框绘制​​
• ​​3. 检测框坐标转换（线矢量）​​
• ​​4. 检测框坐标转换（面矢量）​​
• ​​5. 检测框成果转换​​

1. 在arcgis中新建线矢量

新建线矢量，添加空间参考，例如wgs_1984。

2. 检测框绘制

3. 检测框坐标转换（线矢量）

地理坐标转像素坐标

# -*- coding: utf-8 -*-
from osgeo import ogr
from osgeo import gdal
from osgeo import osr
import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt


def getSRSPair(dataset):
    '''
    获得给定数据的投影参考系和地理参考系
    :param dataset: GDAL地理数据
    :return: 投影参考系和地理参考系
    '''
    prosrs = osr.SpatialReference()
    prosrs.ImportFromWkt(dataset.GetProjection())
    geosrs = prosrs.CloneGeogCS()
    return prosrs, geosrs


def geo2lonlat(dataset, x, y):
    '''
    将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
    :param dataset: GDAL地理数据
    :param x: 投影坐标x
    :param y: 投影坐标y
    :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(prosrs, geosrs)
    coords = ct.TransformPoint(x, y)
    return coords[:2]


def lonlat2geo(dataset, lon, lat):
    '''
    将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）
    :param dataset: GDAL地理数据
    :param lon: 地理坐标lon经度
    :param lat: 地理坐标lat纬度
    :return: 经纬度坐标(lon, lat)对应的投影坐标
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(geosrs, prosrs)
    coords = ct.TransformPoint(lon, lat)
    return coords[:2]


def imagexy2geo(dataset, row, col):
    '''
    根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
    :param dataset: GDAL地理数据
    :param row: 像素的行号
    :param col: 像素的列号
    :return: 行列号(row, col)对应的投影坐标或地理坐标(x, y)
    '''
    trans = dataset.GetGeoTransform()
    px = trans[0] + col * trans[1] + row * trans[2]
    py = trans[3] + col * trans[4] + row * trans[5]
    return px, py


def geo2imagexy(dataset, x, y):
    '''
    根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
    :param dataset: GDAL地理数据
    :param x: 投影或地理坐标x
    :param y: 投影或地理坐标y
    :return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
    '''
    trans = dataset.GetGeoTransform()
    a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
    b = np.array([x - trans[0], y - trans[3]])
    return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解


def main(imgPath, shpPath):
    dataset = gdal.Open(imgPath)
    ds = ogr.Open(shpPath,1)
    if ds is None:
        print('Could not open folder')
    in_lyr = ds.GetLayer()

    lyr_dn = in_lyr.GetLayerDefn()
    cls_index = lyr_dn.GetFieldIndex("Id")
    cls_name_g = lyr_dn.GetFieldDefn(cls_index)
    feature = in_lyr.GetNextFeature()

    fieldName = cls_name_g.GetNameRef()

    finalResult = []
    while feature is not None:
        geom = feature.geometry()
        # ff = feature.GetGeometry(geom)
        # print(geom)
        print('------------')
        result = []
        for i in range(0, geom.GetPointCount()-1):
            pt = np.array(geom.GetPoint(i))
            coords = lonlat2geo(dataset, pt[0], pt[1])
            coords = geo2imagexy(dataset, coords[0], coords[1])
            result.append([coords[0], coords[1]])
            # print(pt, x, y)

        finalResult.append(result)
        feature = in_lyr.GetNextFeature()

    return finalResult


if __name__ == '__main__':
    img_filename = './test/0000000004_image.tif'
    dst_filename = './test/label2.shp'
    finalResult = main(img_filename, dst_filename)
    img = cv.imread(img_filename, cv.IMREAD_LOAD_GDAL)
    finalResult = np.array(finalResult)
    for bbox in finalResult:
        xmin = min(bbox[0][0], bbox[1][0], bbox[2][0], bbox[3][0])
        ymin = min(bbox[0][1], bbox[1][1], bbox[2][1], bbox[3][1])
        xmax = max(bbox[0][0], bbox[1][0], bbox[2][0], bbox[3][0])
        ymax = max(bbox[0][1], bbox[1][1], bbox[2][1], bbox[3][1])
        cv.rectangle(img, (int(xmin), int(ymin)), (int(xmax), int(ymax)), (0, 100, 255), 5)

    plt.imshow(img)
    plt.show()      

4. 检测框坐标转换（面矢量）

大部分代码一样，就是main函数有点区别。

def main(imgPath, shpPath):
    dataset = gdal.Open(imgPath)
    ds = ogr.Open(shpPath,1)
    if ds is None:
        print('Could not open folder')
    in_lyr = ds.GetLayer()

    lyr_dn = in_lyr.GetLayerDefn()
    cls_index = lyr_dn.GetFieldIndex("Id")
    cls_name_g = lyr_dn.GetFieldDefn(cls_index)
    feature = in_lyr.GetNextFeature()

    fieldName = cls_name_g.GetNameRef()

    finalResult = []
    while feature is not None:
        geom = feature.geometry()
        arr = np.array(feature.GetGeometryRef().GetEnvelope())
        print(arr)
        coordsMin = lonlat2geo(dataset, arr[0], arr[3])
        coordsMin = geo2imagexy(dataset, coordsMin[0], coordsMin[1])
        coordsMax = lonlat2geo(dataset, arr[1], arr[2])
        coordsMax = geo2imagexy(dataset, coordsMax[0], coordsMax[1])
        finalResult.append([coordsMin[0], coordsMin[1], coordsMax[0], coordsMax[1]])
        # ff = feature.GetGeometry(geom)
        # print(geom)
        # print('------------')
        # result = []
        # print(geom.GetPointCount())
        # for i in range(0, geom.GetPointCount()-1):
        #     pt = np.array(geom.GetPoint(i))
        #     coords = lonlat2geo(dataset, pt[0], pt[1])
        #     coords = geo2imagexy(dataset, coords[0], coords[1])
        #     result.append([coords[0], coords[1]])
            # print(pt, x, y)

        # finalResult.append(result)
        feature = in_lyr.GetNextFeature()

    return finalResult      

5. 检测框成果转换

# -*- coding: utf-8 -*-
import gdal
from osgeo import ogr, osr


def read_img(filename):
    dataset=gdal.Open(filename)

    im_width = dataset.RasterXSize
    im_height = dataset.RasterYSize

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0,0,im_width,im_height)

    # del dataset
    return im_width,im_height,im_proj,im_geotrans,im_data,dataset


def getSRSPair(dataset):
    '''
    获得给定数据的投影参考系和地理参考系
    :param dataset: GDAL地理数据
    :return: 投影参考系和地理参考系
    '''
    prosrs = osr.SpatialReference()
    prosrs.ImportFromWkt(dataset.GetProjection())
    geosrs = prosrs.CloneGeogCS()
    return prosrs, geosrs


def geo2lonlat(dataset, x, y):
    '''
    将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
    :param dataset: GDAL地理数据
    :param x: 投影坐标x
    :param y: 投影坐标y
    :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(prosrs, geosrs)
    coords = ct.TransformPoint(x, y)
    return coords[:2]


def imagexy2geo(dataset, col, row):
    '''
    根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
    :param dataset: GDAL地理数据
    :param row: 像素的行号
    :param col: 像素的列号
    :return: 行列号(row, col)对应的投影坐标或地理坐标(x, y)
    '''
    trans = dataset.GetGeoTransform()
    print(trans)
    print(row,col)
    px = trans[0] + col * trans[1] + row * trans[2]
    py = trans[3] + col * trans[4] + row * trans[5]
    return px, py


def imagexy2shp(img_path, strVectorFile, bboxes):
    im_width,im_height,im_proj,im_geotrans,im_data,dataset = read_img(img_path)
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO")  # 为了支持中文路径
    gdal.SetConfigOption("SHAPE_ENCODING", "CP936")  # 为了使属性表字段支持中文
    ogr.RegisterAll()
    strDriverName = "ESRI Shapefile"  # 创建数据，这里创建ESRI的shp文件
    oDriver = ogr.GetDriverByName(strDriverName)
    if oDriver == None:
        print("%s 驱动不可用！\n", strDriverName)

    oDS = oDriver.CreateDataSource(strVectorFile)  # 创建数据源
    if oDS == None:
        print("创建文件【%s】失败！", strVectorFile)

    # srs = osr.SpatialReference()  # 创建空间参考
    # srs.ImportFromEPSG(4326)  # 定义地理坐标系WGS1984
    srs = osr.SpatialReference(wkt=dataset.GetProjection())#我在读栅格图的时候增加了输出dataset，这里就可以不用指定投影，实现全自动了，上面两行可以注释了，并且那个proj参数也可以去掉了，你们自己去掉吧
    papszLCO = []
    # 创建图层，创建一个多边形图层,"TestPolygon"->属性表名
    oLayer = oDS.CreateLayer("TestPolygon", srs, ogr.wkbPolygon, papszLCO)
    if oLayer == None:
        print("图层创建失败！\n")

    '''下面添加矢量数据，属性表数据、矢量数据坐标'''
    oFieldID = ogr.FieldDefn("FieldID", ogr.OFTInteger)  # 创建一个叫FieldID的整型属性
    oLayer.CreateField(oFieldID, 1)

    # oFieldName = ogr.FieldDefn("FieldName", ogr.OFTString)  # 创建一个叫FieldName的字符型属性
    # oFieldName.SetWidth(100)  # 定义字符长度为100
    # oLayer.CreateField(oFieldName, 1)

    oDefn = oLayer.GetLayerDefn()  # 定义要素
    # 创建单个面
    for bbox in bboxes:
        oFeatureTriangle = ogr.Feature(oDefn)
        oFeatureTriangle.SetField(0, 0)  # 第一个参数表示第几个字段，第二个参数表示字段的值
        # oFeatureTriangle.SetField(1, "单个面")
        ring = ogr.Geometry(ogr.wkbLinearRing)  # 构建几何类型:线
        ring.AddPoint(bbox[0], bbox[1])  # 添加点01
        ring.AddPoint(bbox[2], bbox[1])  # 添加点02
        ring.AddPoint(bbox[2], bbox[3])  # 添加点03
        ring.AddPoint(bbox[0], bbox[3])  # 添加点04
        yard = ogr.Geometry(ogr.wkbPolygon)  # 构建几何类型:多边形
        yard.AddGeometry(ring)
        yard.CloseRings()

        geomTriangle = ogr.CreateGeometryFromWkt(str(yard))  # 将封闭后的多边形集添加到属性表
        oFeatureTriangle.SetGeometry(geomTriangle)
        oLayer.CreateFeature(oFeatureTriangle)

    oDS.Destroy()
    print("数据集创建完成！\n")


def parse_txt(path):
    f = open(path)
    data = f.readlines()
    anns = []
    for ann in data:
        try:
            ann = ann.split(' ')
            xmin = float(ann[0])
            ymin = float(ann[1])
            xmax = float(ann[2])
            ymax = float(ann[3])
            anns.append([xmin, ymin, xmax, ymax])
        except:
            print('error', ann)
    return anns


if __name__ == "__main__":
    img_filename = './test/0000000004_image.tif'
    dst_filename = './test/label6.shp'
    # imagexy2shp(img_filename, dst_filename)
    txtPath = './test/0000000004_image.txt'
    dataset = gdal.Open(img_filename)
    anns = parse_txt(txtPath)
    annsGEO = []
    for ann in anns:
        xmin, ymin = imagexy2geo(dataset, ann[0], ann[1])
        xmax, ymax = imagexy2geo(dataset, ann[2], ann[3])
        annsGEO.append([xmin, ymin, xmax, ymax])

    # 转换成面矢量
    imagexy2shp(img_filename, dst_filename, annsGEO)