2020-11-17

遥感图像的读取,16位转8位，单通道转三通道，线性拉伸

import os
import numpy as np
from osgeo import gdal

# 定义读取和保存图像，以及仿射变换的类
class GRID:
    def load_image(self, filename):
        image = gdal.Open(filename) #打开图像
 
        img_width = image.RasterXSize
        img_height = image.RasterYSize
        outbandsize = image.RasterCount  #获取图像的宽高及通道数

        img_proj = image.GetProjection()    #获取地图投影信息
        img_geotrans = image.GetGeoTransform() #获取图像仿射矩阵
        img_data = image.ReadAsArray(0, 0, img_width, img_height)   #将数据写成数组，对应栅格矩阵
 
        del image
 
        return outbandsize, img_proj, img_geotrans, img_data
 
    def write_image(self, filename, img_proj, img_geotrans, img_data):
        # 判断栅格数据类型
        if 'uint8' in img_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'uint16' in img_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
 
        # 判断数组维度
        if len(img_data.shape) == 3:
            img_bands, img_height, img_width = img_data.shape
        else:
            img_bands, (img_height, img_width) = 1, img_data.shape
 
        # 创建文件
        driver = gdal.GetDriverByName('GTiff')
        image = driver.Create(filename, img_width, img_height, img_bands, datatype)
 
        image.SetGeoTransform(img_geotrans)
        image.SetProjection(img_proj)
 
        if img_bands == 1:
            image.GetRasterBand(1).WriteArray(img_data)
        else:
            for i in range(img_bands):
                image.GetRasterBand(i+1).WriteArray(img_data[i])
 
        del image # 删除变量,保留数据

    ###仿射矩阵信息有六个参数，描述的是栅格行列号和地理坐标之间的关系####
    '''
    0：图像左上角的X坐标；
    1：图像东西方向分辨率；
    2：旋转角度，如果图像北方朝上，该值为0；
    3：图像左上角的Y坐标；
    4：旋转角度，如果图像北方朝上，该值为0；
    5：图像南北方向分辨率；
    '''
    def CoordTransf(self,Xpixel,Ypixel,GeoTransform):
        XGeo = GeoTransform[0]+GeoTransform[1]*Xpixel+Ypixel*GeoTransform[2]
        YGeo = GeoTransform[3]+GeoTransform[4]*Xpixel+Ypixel*GeoTransform[5]
        return XGeo,YGeo
 
# 图像对比度增强之截断线性拉伸
class PercentLinear:
	# 单通道图像2%线性拉伸
    def gray_process(self, gray, maxout = 255, minout = 0):
        high_value = np.percentile(gray, 98)#取得98%直方图处对应灰度
        low_value = np.percentile(gray, 2)#取得2%直方图处对应灰度
        truncated_gray = np.clip(gray, a_min=low_value, a_max=high_value) 
        processed_gray = ((truncated_gray - low_value)/(high_value - low_value)) * (maxout - minout)#线性拉伸
        return np.uint8(processed_gray)
	
	# 三通道图像2%线性拉伸
    def rgb_process(self, image, max_out=255, min_out=0):
        b, g, r = cv2.split(image)#分开三个波段
        r_p = gray_process(r)
        g_p = gray_process(g)
        b_p = gray_process(b)
        result = cv2.merge((b_p, g_p, r_p))#合并处理后的三个波段
        return result

1. Gdal处理遥感图像

#读取图像
path = 'read.tif'
bandsize, proj, geotrans, data = GRID().load_image(path)

#保存图像
GRID().write_image('save.tif', proj, geotrans, data)

#仿射变换
img_height, img_width = data.shape
GRID().CoordTransf(img_height,img_width,geotrans)

2. 将16位图像转为8位

data = np.uint8(data)

3. 将单通道图像转为三通道

data = cv2.merge((data, data, data))

4. 对图像进行线性拉伸

data = PercentLinear().gray_process(data)
data = PercentLinear().rgb_process(data)