欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

如何利用Python对遥感影像进行显示

程序员文章站 2022-03-22 15:00:34
遥感影像一般具有多个波段,比较常见的影像一般是4波段多光谱影像,比如高分一号、高分二号、资源三号等。这些影像数据一般体量较大,有的几百兆,有的多达几十G,格式一般是16位无符号整型,一般看图软件无法打开显示,需要ArcGIS、ENVI等专业的软件进行查看,有时候很不方便。这篇博客就简单的介绍一下,如何利用Python对遥感影像进行显示,需要用到的库为GDAL和Opencv。正文Python中,一般的绘图库都无法处理或显示遥感影像,例如Matplotlib、Opencv、Scipy等,很多连遥...

遥感影像一般具有多个波段,比较常见的影像一般是4波段多光谱影像,比如高分一号、高分二号、资源三号等。

这些影像数据一般体量较大,有的几百兆,有的多达几十G,格式一般是16位无符号整型,一般看图软件无法打开显示,需要ArcGIS、ENVI等专业的软件进行查看,有时候很不方便。

这篇博客就简单的介绍一下,如何利用Python对遥感影像进行显示,需要用到的库为GDAL和Opencv。

正文

Python中,一般的绘图库都无法处理或显示遥感影像,例如Matplotlib、Opencv、Scipy等,很多连遥感影像都读取不了。

基于此,本代码的主要解决思路:(1)利用GDAL读取遥感影像,获取影像数据;(2)对影像数据进行简单的拉伸变换处理,处理成8位整型;(3)将处理后的数据显示出来。附带增加了一个Gamma变换函数,可以对数据进行拉伸处理。

测试数据:高分一号8米4波段多光谱数据

这里,遥感影像的读取函数,来源于之前写的博客 如何使用Python中的GDAL库对遥感影像进行读取和存储。

# -*- coding: utf-8 -*-
import numpy as np
from osgeo import gdal
import cv2

class IMAGE:
    #读图像文件
    def read_img(self, filename):
        dataset = gdal.Open(filename) #打开文件
        im_width = dataset.RasterXSize  #栅格矩阵的列数
        im_height = dataset.RasterYSize  #栅格矩阵的行数
        im_bands = dataset.RasterCount   #波段数
        im_geotrans = dataset.GetGeoTransform()  #仿射矩阵,左上角像素的大地坐标和像素分辨率
        im_proj = dataset.GetProjection() #地图投影信息,字符串表示
        im_data = dataset.ReadAsArray(0,0,im_width,im_height)
 
        del dataset 
 
        return im_bands, im_data
 
    #写GeoTiff文件
    def write_img(self, filename, im_proj, im_geotrans, im_data):
 
        #判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
 
        #判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
           im_bands, (im_height, im_width) = 1, im_data.shape
 
        #创建文件
        driver = gdal.GetDriverByName("GTiff") 
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
 
        dataset.SetGeoTransform(im_geotrans)       #写入仿射变换参数
        dataset.SetProjection(im_proj)          #写入投影
 
        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i+1).WriteArray(im_data[i])
 
        del dataset

## 影像拉伸 ##
def img_processing(im_band,img_data):
    if im_band == 1:
        data_jpg = np.zeros((img_data.shape[0],img_data.shape[1]),dtype='uint8')
        im_max = np.amax(img_data)
        im_min = np.amin(img_data)
        for m in range(0, img_data.shape[0]):
            for n in range(0, img_data.shape[1]):
                data_jpg[m,n] = 255/(im_max-im_min)*(img_data[m,n]-im_min)
    else:
        data_jpg = np.zeros((img_data.shape[1],img_data.shape[2],3),dtype='uint8')
        for i in range(3):
            im_max = np.amax(img_data[i,:,:])
            im_min = np.amin(img_data[i,:,:])
            for m in range(0, img_data.shape[1]):
                for n in range(0, img_data.shape[2]):
                    data_jpg[m,n,i] = 255/(im_max-im_min)*(img_data[i,m,n]-im_min)
    return data_jpg

## GAMMA 变换 ##
def gamma_trans(img, gamma):
    gamma_table = [np.power(x/255.0,gamma)*255.0 for x in range(256)]
    gamma_table = np.round(np.array(gamma_table)).astype(np.uint8)
    return cv2.LUT(img,gamma_table)

if __name__ == '__main__':
    image_grid = IMAGE()
    image_name = 'GF1.tif'
    im_band,img_data = image_grid.read_img(image_name)
    data_jpg = img_processing(im_band,img_data)
    data_jpg_corrted = gamma_trans(data_jpg, 0.5)
    ## 显示 ##
    cv2.imshow('GF1',data_jpg)
    cv2.imshow('Gamma_image',data_jpg_corrted)
    cv2.waitKey()

ArcGIS下显示

拉伸参数设置:如何利用Python对遥感影像进行显示

如何利用Python对遥感影像进行显示

经过Python拉伸处理后的图

如何利用Python对遥感影像进行显示

经过Gamma变换(系数为0.5)后的图

如何利用Python对遥感影像进行显示

总结

可以看出,处理后的效果和ArcGIS中的显示效果还是有较大差距,不过图像内容还是能够看出来。其他数据由于没有,所以暂时还没有测试,可能会存在问题,不是可能,严谨点,肯定存在问题。后续根据需要,再慢慢解决吧!

本文地址:https://blog.csdn.net/u012569919/article/details/109864492