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

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

文章目录

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

1. 在arcgis中新建线矢量

2. 检测框绘制

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

地理坐标转像素坐标

# -*- coding: utf-8 -*-from osgeo import ogrfrom osgeo import gdalfrom osgeo import osrimport numpy as npimport cv2 as cvimport matplotlib.pyplot as pltdef getSRSPair(dataset): ''' 获得给定数据的投影参考系和地理参考系 :param dataset: GDAL地理数据 :return: 投影参考系和地理参考系 ''' prosrs = osr.SpatialReference() prosrs.ImportFromWkt(dataset.GetProjection()) geosrs = prosrs.CloneGeogCS() return prosrs, geosrsdef 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, pydef 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 finalResultif __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 gdalfrom osgeo import ogr, osrdef 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,datasetdef getSRSPair(dataset): ''' 获得给定数据的投影参考系和地理参考系 :param dataset: GDAL地理数据 :return: 投影参考系和地理参考系 ''' prosrs = osr.SpatialReference() prosrs.ImportFromWkt(dataset.GetProjection()) geosrs = prosrs.CloneGeogCS() return prosrs, geosrsdef 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, pydef 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 annsif __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)

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。