基于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)