Python 遥感图像读写 python处理遥感影像

2023-06-27 21:32| 来源: 网络整理| 查看: 265

遥感影像一般具有多个波段，比较常见的影像一般是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] = float(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] = float(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处理遥感影像_opencv

Python 遥感图像读写 python处理遥感影像_numpy_02

经过Python拉伸处理后的图

Python 遥感图像读写 python处理遥感影像_opencv_03

经过Gamma变换（系数为0.5）后的图

Python 遥感图像读写 python处理遥感影像_numpy_04

总结

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

【本文地址】

公司简介

联系我们