（Python）使用Gdal+opencv读遥感影像加简单处理

简介

使用Python API的gdal和opencv一般遇到最大的问题就是，遥感影像一般为16位和opencv只处理8位，中间的转换有很多的写的非常复杂又不完全对的，本文提供一个可用的例子，用于读取影像截取ROI，这样就可以很方便的使用opencv算法进行图像处理了，以此作为自己的记录和初学者参考。

使用库的相关版本
python 3.9
gdal 3.2.3
opencv 4.5（可使用PLI或者matplot）

代码

from osgeo import gdal
import numpy as np
import cv2
from time import *
def readTIFF(tifpath, bandnum):
    """
    Use GDAL to read data and transform them into arrays.
    :param tifpath:tif文件的路径
    :param bandnum:需要读取的波段
    :return:该波段的数据，narray格式。len(narray)是行数，len(narray[0])列数
    """
    image = gdal.Open(tifpath)  # 打开影像
    if image == None:
        print(tifpath + "该tif不能打开!")
        return
    im_width = image.RasterXSize  # 栅格矩阵的列数，经度方向，cols(好多人搞反)
    im_height = image.RasterYSize  # 栅格矩阵的行数，纬度方向，rows
    im_bands = image.RasterCount  # 波段数
    im_proj = image.GetProjection()  # 获取投影信息坐标系
    im_geotrans = image.GetGeoTransform()  # 仿射矩阵，具体参考GetGeoTransform
    print('tif数据:{}个行，{}个列，{}层波段, 取出第{}层.'.format(im_width, im_height, im_bands, bandnum))
    im_data = image.ReadAsArray(0, 0,  im_width, im_height)
    del image  # 减少冗余
    return im_data,im_proj, im_geotrans
def normalization(data):
    _range = np.max(data) - np.min(data)
    return (data - np.min(data)) / _range
def Tiff16to8bit(img_16):
    if (np.max(img_16) - np.min(img_16) != 0):
        # img_nrm = (img_16 - np.min(img_16)) / (np.max(img_16) - np.min(img_16)) #计算灰度范围,归一化
        img_nrm = normalization(img_16)
        img_8 = np.uint8(255 * img_nrm)
    return img_8
begin_time = time()
#读入图片
im_data,im_proj, im_geotrans = readTIFF("D:\\SongshanZY3\\ortho\\BWDSC.tif",0)
#取小块的左上角坐标
x_offset,y_offset = 5000,6000
#小块大小
block_size_x = 1000
block_size_y = 2000
#注意矩阵访问和直接的坐标系是相反的，这个在C++的at<double>(y,x)也有体现
block_im = im_data[y_offset:y_offset + block_size_y,x_offset:x_offset + block_size_x]
block_im = Tiff16to8bit(block_im)
end_time = time()
run_time = end_time-begin_time
print ('该程序运行时间：',run_time,'s') #该循环程序运行时间： 1.4201874732
#展示并存储
cv2.imwrite("block_B.png",block_im)
cv2.namedWindow("block",0)#创建可以缩放大小的窗口
cv2.imshow("block",block_im)

#进行一些opencv的测试使用
#canny边缘提取
block_im = cv2.GaussianBlur(block_im, (3, 3), 0)
canny = cv2.Canny(block_im, 30, 120)
cv2.namedWindow("Canny",0)#创建可以缩放大小的窗口
cv2.imshow('Canny', canny)

#hough直线检测
hough_img = block_im.copy()
lines = cv2.HoughLinesP(canny,1,np.pi/180,30,minLineLength=60,maxLineGap=10)
lines1 = lines[:,0,:]#提取为二维
for x1,y1,x2,y2 in lines1[:]:
    cv2.line(hough_img,(x1,y1),(x2,y2),(255,0,0),1)
cv2.namedWindow("Hough",0)#创建可以缩放大小的窗口
cv2.imshow("Hough", hough_img)

#otsu阈值分割
ret2, otsu = cv2.threshold(block_im, 0 , 255, cv2.THRESH_OTSU)
cv2.namedWindow("Otsu",0)#创建可以缩放大小的窗口
cv2.imshow("Otsu", otsu)

cv2.waitKey()

运行结果

可以看到30000*20000的图在笔者的thinkpadX1上时间基本1s左右，速度尚可；python的gdal比C++写起来方便多了，但是接口没那么全。
而且可以看到Opencv可以快乐的使用起来了，希望为初学者带来一些参考，后续继续补充。如需交流，个人微信WHUwsd1995.