得到黑白结果图后
由于某些原因需要将结果图转为矢量图进行后处理
这是由于大多数软件比如arcgis对矢量shp文件处理较为方便
而对栅格图片的处理有很多不便
因此我们需要将我们的黑白栅格图片转为矢量shp文件
这里主要使用arcgis软件实现
若直接将深度学习的结果图片导入会发现它没有地理坐标！！
与原图不重合
那么矢量化的面状则无法合并
这是由于我的深度学习代码储存图片使用的是cv2库
这个库不会保留原图的地理坐标
因此需要使用gdal处理遥感图像的库来实现我们的需求
所以我们首先需要将原图的地理坐标赋给我们生成的png深度学习结果图
代码实现——将有坐标的遥感图像的地理坐标赋给不带有坐标的png图
在这里希望大家还能记得我们深度学习的样本需要裁剪成512*512大小
但是通过cv2这个库实现的图片不再保留原图的地理坐标
这部分代码可以查看我的这篇博文
但是我现在想在裁剪的同时保留对应的地理坐标和投影信息
那该怎么修改代码呢
这就需要使用处理遥感影像的gdal库了
关于gdal库的一些接口可以查看这篇博文
首先我的代码可以实现对包含投影信息和地理坐标的tif图进行裁剪处理
同时对包含地理信息的png图或者不包含地理信息的png图也可以进行裁剪处理(这部分是arcgis软件做出来的）
最后我还是添加了筛选删除掉全黑背景样本的代码
可以直接把没有道路的图片删除掉
话不多说
直接上我实现这个功能的代码
# -*- coding: utf-8 -*-
"""
Created on Wed July 19 11:44:46 2022
@author:Laney_Midory
csdn:Laney_Midory
"""
import cv2
import os
from osgeo import gdal
# Cutting the input image to h*w blocks
#heightCutNum = 2
#widthCutNum = 2
# The folder path of input and output
#inPath = "C:/Users/Administrator/Desktop/train-water-png/"
inPath = "C:/Users/Administrator/Desktop/train-road/"
outPath = "C:/Users/Administrator/Desktop/out/"
#outPath = "C:/Users/Administrator/Desktop/t/"
for f in os.listdir(inPath):
path = inPath + f.strip()
print(path)
#print(outPath + str(1) + str(2) + f.strip())
in_ds = gdal.Open(path)
# 读取原图中的每个波段，通道数从1开始，默认前三波段
in_band1 = in_ds.GetRasterBand(1)
in_band2 = in_ds.GetRasterBand(2)
in_band3 = in_ds.GetRasterBand(3)
# 获取原图的原点坐标信息
ori_transform = in_ds.GetGeoTransform()
top_left_x = ori_transform[0] # 左上角x坐标
top_left_y = ori_transform[3] # 左上角y坐标
w_e_pixel_resolution = ori_transform[1] # 东西方向像素分辨率
n_s_pixel_resolution = ori_transform[5] # 南北方向像素分辨率
# The size of each input image
height = in_ds.RasterYSize # 栅格矩阵的行数
width = in_ds.RasterXSize # 栅格矩阵的列数
# The size of block that you want to cut
#heightBlock = int(height / heightCutNum)
#widthBlock = int(width / widthCutNum)
offset_x = 0
offset_y = 0
block_xsize = 512
block_ysize = 512
num_width = int(width / block_ysize)
num_height = int(height / block_xsize)
count = 0
for i in range(num_width):
for j in range(num_height):
count = count + 1
offset_x1 = offset_x + block_xsize * i
offset_y1 = offset_y + block_ysize * j
out_band1 = in_band1.ReadAsArray(offset_x1, offset_y1, block_xsize, block_ysize)
out_band2 = in_band2.ReadAsArray(offset_x1, offset_y1, block_xsize, block_ysize)
out_band3 = in_band3.ReadAsArray(offset_x1, offset_y1, block_xsize, block_ysize)
gtif_driver = gdal.GetDriverByName("GTiff") # 数据类型必须有，因为要计算需要多大内存空间，但是这儿是只有GTiff吗？
filename = outPath + str(count) + f.strip() # 文件名称
out_ds = gtif_driver.Create(filename, block_xsize, block_ysize, 3, in_band1.DataType) # 数据格式遵循原始图像
top_left_x1 = top_left_x + offset_x1 * w_e_pixel_resolution
top_left_y1 = top_left_y + offset_y1 * n_s_pixel_resolution
dst_transform = (top_left_x1, ori_transform[1], ori_transform[2], top_left_y1, ori_transform[4], ori_transform[5])
out_ds.SetGeoTransform(dst_transform)
# 设置SRS属性（投影信息）
out_ds.SetProjection(in_ds.GetProjection())
# 写入目标文件（如果波段数有更改，这儿也需要修改）
out_ds.GetRasterBand(1).WriteArray(out_band1)
out_ds.GetRasterBand(2).WriteArray(out_band2)
out_ds.GetRasterBand(3).WriteArray(out_band3)
# 将缓存写入磁盘，直接保存到了程序所在的文件夹
# out_ds.FlushCache()
print(f"FlushCache succeed{count}")
del out_ds
这部分为裁剪并且添加上地理坐标
下面实现删除掉黑色背景的图片
for f in os.listdir(outPath):
path = outPath + f.strip()
if not os.path.exists(path):
continue;
img = cv2.imread(path,0)
if cv2.countNonZero(img) == 0:
print("Image is black")
print(path)
path2=f.strip().split("_")
print(outPath +path2[0] + "_sat.tif")
os.remove(path)
os.remove(outPath +path2[0] + "_sat.tif")
print("finish!")
这两部分拼接起来运行速度其实还快于cv2的处理速度
我使用gdal库以及cv2库对一样的图片进行处理
最后都得到相同数量的裁剪图片
说明代码效果没问题
而速度gdal还快于cv2
但这里还有个问题
由于我的png图片不全是RGB图像
部分图片位深度是8位，部分是24位
我上面的代码处理的是24位的图片
否则遇到8位图片会报错
Traceback (most recent call last):
File “c:/Users/Administrator/Desktop/DeepGlobe-Road-Extraction-link34/im8to24.py”, line 30, in
print(in_ds.getbands())
AttributeError: ‘Dataset’ object has no attribute ‘getbands’
于是对代码进行进一步修改
在这里我print出来，看看是哪一个需要修改
运行结果是
这样可以看出确实是图片的问题
需要进行进一步处理
首先我考虑了直接把8位深度转为24位深度
但是gdal库没有直接的函数能实现（可能是我没找到，希望有知道的朋友分享一下）
我搜索了之后发现这篇博文使用了这个方法
但是我担心这样一转之后就不包含地理信息了
所以最好不进行这个处理
于是我采取了第二种方法
这里就需要在源代码的基础上加上对于深度为1的情况
代码如下
if bands_num == 1:
#in_ds.GetRasterBand(1).WriteArray(img_array)
in_band1 = in_ds.GetRasterBand(1)
else:
in_band1 = in_ds.GetRasterBand(1)
in_band2 = in_ds.GetRasterBand(2)
in_band3 = in_ds.GetRasterBand(3)
在源代码的基础上
需要以类似的方法修改三处
最好别调换顺序
我就是因为调换了顺序出现了点问题
由此就可以顺利解决所有情况下的裁剪并保留地理坐标的问题啦！
arcgis软件实现——栅格数据矢量化
线要素较多的情况下可以查看这篇文章
但我们的图这种情况可以使用这个方法
首先是导入我们有地理位置的图片
点击栅格转面可以直接按照value转为面状
但这样结果是整个面
打开编辑器
可以看见很多个面
我们只需要把我们不需要的面删除即可
这样就可以得到我们需要的面啦
将裁剪之后的原图坐标赋给对应的深度学习结果图
由于Dinknet深度学习保存和处理图片使用的是cv2库
因此处理之后获取的深度学习结果图是不带有地理坐标信息的
经过考虑之后
对深度学习代码将cv2修改为gdal库有些困难
因此这里还是对处理之后的结果图进行后处理
即再进行一步将裁剪之后的原图坐标赋给对应的深度学习结果图
接口跟第一步是相同的
但这里是直接把一张图片的地理信息赋给另外一张图片
因此需要读取两张图片
一张是有地理信息的
一张是不含有地理信息的
代码如下
# -*- coding: utf-8 -*-
"""
Created on Wed July 20 12:01:20 2022
@author:Laney_Midory
csdn:Laney_Midory
"""
import os
from osgeo import osr, gdal
def assign_spatial_reference_byfile(src_path, dst_path):
src_ds = gdal.Open(src_path, gdal.GA_ReadOnly)
print(type(src_ds))
sr = osr.SpatialReference()
sr.ImportFromWkt(src_ds.GetProjectionRef())
geoTransform = src_ds.GetGeoTransform()
dst_ds = gdal.Open(dst_path, gdal.GA_Update)
print(type(dst_ds))
dst_ds.SetProjection(sr.ExportToWkt())
dst_ds.SetGeoTransform(geoTransform)
dst_ds = None
src_ds = None
'''
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326) #地理坐标投影
ds = gdal.Open(dst_path, gdal.GA_Update)
ds.SetProjection(sr.ExportToWkt())
ds.SetGeoTransform([0,1,0,0,0,1]) #地理6参数
ds = None
'''
def def_geoCoordSys(read_path, img_transf, img_proj):
array_dataset = gdal.Open(read_path)
img_array = array_dataset.ReadAsArray(
0, 0, array_dataset.RasterXSize, array_dataset.RasterYSize)
if 'int8' in img_array.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in img_array.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(img_array.shape) == 3:
img_bands, im_height, im_width = img_array.shape
else:
img_bands, (im_height, im_width) = 1, img_array.shape
filename = read_path[:-4] + '_proj' + '.tif'
print(filename)
driver = gdal.GetDriverByName("GTiff") # 创建文件驱动
dataset = driver.Create(
filename, im_width, im_height, img_bands, datatype)
dataset.SetGeoTransform(img_transf) # 写入仿射变换参数
dataset.SetProjection(img_proj) # 写入投影
# 写入影像数据
if img_bands == 1:
dataset.GetRasterBand(1).WriteArray(img_array)
else:
for i in range(img_bands):
dataset.GetRasterBand(i + 1).WriteArray(img_array[i])
print(read_path, 'geoCoordSys get!')
src_path = 'C:/Users/Administrator/Desktop/2/1ad04_sat.tif'
dst_path = 'C:/Users/Administrator/Desktop/2/00ad04_mask.png'
#assign_spatial_reference_byfile(src_path, dst_path)
dataset = gdal.Open(src_path) # 打开文件
img_pos_transf = dataset.GetGeoTransform() # 仿射矩阵
img_pos_proj = dataset.GetProjection() # 地图投影信息
def_geoCoordSys(dst_path,img_pos_transf, img_pos_proj)
print("finish!")
其中src_path是包含有地理信息的图片
而dst_path是需要修改的目标图片
但这样一张张处理有些复杂
因此把上面的代码改成批量处理
并且设置一个目标文件夹存储我们处理好的文件
因此先对我们的函数加一个存储文件夹
def def_geoCoordSys(read_path, img_transf, img_proj,out_path):
对存储路径进行修改
#filename = read_path[:-4] + '_proj' + '.tif'
filename = out_path+'\\'+file[:-4] + '_proj' + '.tif'
print(filename)
主函数再进行一个修改
change_path = r"C:\Users\Administrator\Desktop\result" # 无地理信息的文件夹
unchange_path = r"C:\Users\Administrator\Desktop\real" # 有地理信息的文件夹
out_path=r"C:\Users\Administrator\Desktop\proj" # 存储文件夹
change_list = os.listdir(change_path)
# print(file_list)
unchange_list = os.listdir(unchange_path)
for file in change_list:
dst_path = os.path.join(change_path, file)#无地理坐标的图片
src_path = os.path.join(unchange_path, file) #有地理坐标的图片
#assign_spatial_reference_byfile(src_path, dst_path)
if not os.path.isfile(unchange_path+'\\'+file):
print("can't find"+unchange_path+'\\'+file)
dataset = gdal.Open(src_path) # 打开文件
img_pos_transf = dataset.GetGeoTransform() # 仿射矩阵
img_pos_proj = dataset.GetProjection() # 地图投影信息
def_geoCoordSys(dst_path,img_pos_transf, img_pos_proj,out_path)
print("finish!")
这样就可以实现批处理啦
处理之后就获取到了后缀为proj的带有地理信息的深度学习结果图
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
