如何利用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拉伸处理后的图
经过Gamma变换(系数为0.5)后的图
总结
可以看出,处理后的效果和ArcGIS中的显示效果还是有较大差距,不过图像内容还是能够看出来。其他数据由于没有,所以暂时还没有测试,可能会存在问题,不是可能,严谨点,肯定存在问题。后续根据需要,再慢慢解决吧!
本文地址:https://blog.csdn.net/u012569919/article/details/109864492
下一篇: 当前云计算行业的落地着力点到底在哪?