前言
遥感图像往往尺寸较大,无法用默认的图像浏览器加载。
GDAL是空间数据处理的开源包,支持多种数据格式的读写。
遥感图像是一种带大地坐标的栅格数据,因此,可以借用GDAL对遥感图像进行读写,本文就来记录一些相关操作。
GDAL的安装和引入
gdal可通过荧光动力学实验室(Laboratory for Fluorescence Dynamics)提供的镜像网站下载安装:
网站链接:https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
有些老版本gdal的引入方式是直接import:
import gdal
新版本的gdal引入方式如下:
from osgeo import gdal
行列数和波段数
下面的示例读取了一张tif遥感图片,输出该栅格数据的行列数和波段数:
from osgeo import gdal
data = gdal.Open("xdu.tif")
rows = data.RasterYSize
cols = data.RasterXSize
bands = data.RasterCount
print(f"rows:{rows}")print(f"cols:{cols}")print(f"bands:{bands}")
输出:
rows:37787
cols:36805
bands:4
坐标转换参数
GetGeoTransform()
方法返回栅格数据的坐标转换参数,即行列坐标与空间坐标的转换参数,示例:
from osgeo import gdal
data = gdal.Open("xdu.tif")
geotrans = data.GetGeoTransform()print(geotrans)
输出:
(298735.10954000003, 0.057460000000000004, 0.0, 3779222.4793800004, 0.0, -0.057460000000000004)
输出值为一个元组,6个元素的含义如下:
- 298735.10954000003:左上角像元x坐标
- 0.057460000000000004:x方向比例尺(像元宽度)
- 0.0:x方向旋转角度
- 3779222.4793800004:左上角像元y坐标
- 0.0:y方向旋转角度
- -0.057460000000000004:y方向比例尺(像元高度)
若影像不含地理坐标信息,默认返回值是:(0,1,0,0,0,1)
空间参照系统信息
·GetProjection()方法返回栅格数据的坐标转换参数,示例:
from osgeo import gdal
data = gdal.Open("xdu.tif")
proj = data.GetProjection()print(proj)
输出:
PROJCS[“WGS 84 / UTM zone 49N”,GEOGCS[“WGS 84”,DATUM[“WGS_1984”,SPHEROID[“WGS 84”,6378137,298.257223563,AUTHORITY[“EPSG”,“7030”]],AUTHORITY[“EPSG”,“6326”]],PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,“8901”]],UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,“9122”]],AUTHORITY[“EPSG”,“4326”]],PROJECTION[“Transverse_Mercator”],PARAMETER[“latitude_of_origin”,0],PARAMETER[“central_meridian”,111],PARAMETER[“scale_factor”,0.9996],PARAMETER[“false_easting”,500000],PARAMETER[“false_northing”,0],UNIT[“metre”,1,AUTHORITY[“EPSG”,“9001”]],AXIS[“Easting”,EAST],AXIS[“Northing”,NORTH],AUTHORITY[“EPSG”,“32649”]]
栅格数据转数组
ReadAsArray()
方法可实现将栅格数据转换成数组(Array)形式,以便后续处理,示例:
from osgeo import gdal
data = gdal.Open("xdu.tif")
data_array = data.ReadAsArray()print(data_array.shape)
band1 = data.GetRasterBand(1)# 获取第一个波段的数据
band1_array = band1.ReadAsArray()print(band1_array.shape)
输出:
(4, 37787, 36805)
(37787, 36805)
对于单波段栅格数据,ReadAsArray()函数返回(rows, columns)
对于多波段栅格数据,ReadAsArray()函数返回(bands, rows, columns)
按块读取栅格
ReadAsArray
同样支持按块读取栅格信息,即读取部分区域图像信息,示例:
from osgeo import gdal
data = gdal.Open("xdu.tif")
data_array = data.ReadAsArray(xoff=15000, yoff=15000, xsize=5, ysize=5)print(data_array)
输出:
[[[ 56 73 64 57 67]
[ 40 48 60 57 64]
[ 35 45 76 65 62]
[ 42 73 93 74 67]
[ 58 80 92 76 74]]
[[ 71 89 81 76 86]
[ 55 63 77 75 84]
[ 51 61 92 83 81]
[ 58 89 110 93 86]
[ 74 96 109 95 94]]
[[ 34 51 41 35 43]
[ 21 28 40 37 43]
[ 19 28 59 46 41]
[ 27 58 76 55 44]
[ 42 64 72 54 49]]
[[255 255 255 255 255]
[255 255 255 255 255]
[255 255 255 255 255]
[255 255 255 255 255]
[255 255 255 255 255]]]
输出四组5x5的矩阵,表示各波段数据。
其中,该函数具体的参数含义如下:
- xoff,yoff:想要读取的部分原点位置在整张图像中距离全图原点的位置
- xsize和ysize指定要读取部分图像的矩形大小
实现大图显示
有些遥感影像地图通常较大,用微软默认的图片查看器无法打开显示。通是借助QGIS、ENVI这类专业软件进行查看,这类软件的显示逻辑基本上是“分层动态加载”,即全局显示时显示缩略图,放大显示时,重新加载局部的精细图,不过存在的问题是浏览不流畅,每次拖动或缩放时,图片均需要消耗时间来进行重新加载。于是思考,有无其它解决方式?
方案一:拉伸变换
图像无法加载的主要原因是加载图像时,需要将图像的每个像素点信息加载进内存,如果将每个像素点所需内存体积减小,就可能能够直接进行加载查看。
这篇博文[3]采用了对图像进行拉伸变化的思路,对图像的每个像素点进行拉伸变换,处理成8位整型。不过经我实测发现,对于大型遥感图像所起到效果有限,并且十分耗时。
方案二:瓦片显示
瓦片是一个遥感术语,是指将一定范围内的地图按照一定的尺寸和格式,切成若干行和列的正方形栅格图片。整幅图显示不了,那就切分成多个瓦片进行分块显示,再进行组装,可以有效减小资源依赖。
这篇博文[4]采用了该方案进行图像显示。经实测,该方案能够有效解决遥感大图显示问题,并且拖动浏览较为流畅,但在显示之前需要耗费一定时间来切分瓦片。下面是瓦片显示实现的核心代码。
- 切分瓦片 第一步是设定瓦片尺寸为1000x1000,然后根据图像的尺寸来进行切片,主要方式是通过
ReadAsArray(w_range_start, h_range_start, w_range, h_range)
来读取区域栅格图像数据,然后保存进字典。
classTiles:def__init__(self, dataset, size=1000, image_type=None):"""
初始化影像瓦片
:param dataset: 影像源数据集
:param size:
"""# 瓦片大小,默认1000
self.size =1000
self.image_type = image_type
# 实际的瓦片字典
self.tiles_dict_source ={}# 显示的瓦片
self.tiles_dict_show ={}
self.dataset = dataset
self.bands = self.dataset.RasterCount
self.width = self.dataset.RasterXSize
self.height = self.dataset.RasterYSize
self.size = size
self.max_value =-99999
self.min_value =99999# 横向瓦片个数
self.w_t = math.ceil(self.width / self.size)# 纵向瓦片个数(图片宽度/瓦片大小(1000))(向上取整数)
self.h_t = math.ceil(self.height / self.size)# 初始化图层极值
self.native_extremum_array = np.zeros([self.bands,2])
self.extremum_array = np.zeros([self.bands,2])
self.static ={}# 获取各个图层的最大和最小值for band inrange(self.bands):
self.native_extremum_array[band]=[-999999999,999999999]
self.extremum_array[band]=[-999999999,999999999]# 初始化瓦片
self.__init_tiles()def__init_tiles(self):"""
瓦片切片
:return:
"""# 遍历纵向瓦片个数for h inrange(self.h_t):
h_range_start = h * self.size # h下表索引,h_range_start初始范围
h_range = self.size
# 最终越界处理if h_range_start + h_range > self.height:
h_range = self.height - h_range_start
# 遍历横向瓦片个数for w inrange(self.w_t):
w_range_start = w * self.size
w_range = self.size
# 最终越界处理if w_range_start + w_range > self.width:
w_range = self.width - w_range_start
# 读取栅格范围数据,保存进字典
tiles = self.dataset.ReadAsArray(w_range_start, h_range_start, w_range, h_range)
self.tiles_dict_source[(h, w)]= tiles
System.signal.signal_progress.emit("正在切分瓦片:",(h +1)*100/ self.h_t)
QApplication.processEvents()
- 影像分析 影像分析这步主要是用来统计更新中的影像极值,以便后续拉伸处理。
defimage_extremum(self):
length =len(self.tiles_dict_source)
cur =0for key, value in self.tiles_dict_source.items():# 统计极值for band inrange(self.bands):
layer = value[band]
max_new = np.max(layer)
min_new = np.min(layer)
max_old = self.native_extremum_array[band][0]
min_old = self.native_extremum_array[band][1]if max_new > max_old:
max_old = max_new
if min_new < min_old:
min_old = min_new
self.native_extremum_array[band]=[max_old, min_old]
cur = cur +1
System.signal.signal_progress.emit("影像分析:",(cur / length)*100)
QApplication.processEvents()
- 初始化影像显示 此步主要是对影像进行拉伸变换,将影像值拉伸到0-255之间,然后将处理后的数据保存进一个新的字典
tiles_dict_show
defstretch_extremum(self, array):"""
影像拉伸-将影像值拉伸到0-255之间
:param array:
:return:
"""
show_tiles = np.empty(array.shape, dtype='uint8')for band inrange(self.bands):# 影像极值
max_value = self.extremum_array[band][0]
min_value = self.extremum_array[band][1]
coefficient =(max_value - min_value)/255if self.bands ==1:
coefficient_array = array[:,:]/ coefficient
else:
coefficient_array = array[band,:,:]/ coefficient
coefficient_array = coefficient_array - min_value / coefficient
show_layer = coefficient_array.astype('uint8')if self.bands ==1:
show_tiles = show_layer
else:
show_tiles[band]= show_layer
if self.bands ==1:return show_tiles
else:return np.rollaxis(show_tiles,0,3)definit_show_tiles(self, native=True, max_pix=255, min_pix=0):"""
初始化用于显示的影像瓦片
:param min_pix:
:param max_pix:
:param native:
:return:
"""
self.__set_extremum(native, max_pix, min_pix)
length =len(self.tiles_dict_source)# 瓦片字典的长度
cur =0for key, value in self.tiles_dict_source.items():# 显示光学影像# 将需要显示的保存到tiles_dict_show字典中if System.m_image_type_optical == self.image_type:
self.tiles_dict_show[key]= self.stretch_extremum(value)# 显示雷达影像elif System.m_image_type_radar == self.image_type:
self.tiles_dict_show[key]= self.radar_show(value,255,0)# 更新进度条信息
cur = cur +1
System.signal.signal_progress.emit("初始化影像显示:",(cur / length)*100)
QApplication.processEvents()
- 创建瓦片并显示 最后这步是根据每一个瓦片的数据,单独创建一个pixmap,从而能够使其在
QGraphicsView
进行显示。
defcreate_pixmap(self, image_array):"""
根据给定的数组生成pixmap
:param image_array:
:return:
"""
shape = image_array.shape
height =0
width =0
bytesPerComponent =1iflen(shape)==2:
height, width = shape
eliflen(shape)==3:
height, width, bytesPerComponent = shape
bytesPerLine = bytesPerComponent * width
if bytesPerComponent ==1:# 显示灰度图像
QImg = QImage(image_array.tobytes(), width, height, bytesPerLine, QImage.Format_Grayscale8)elif bytesPerComponent ==3:# 显示GGB图像
QImg = QImage(image_array.tobytes(), width, height, bytesPerLine, QImage.Format_RGB888)else:# 如果不符合就不显示returnreturn QPixmap.fromImage(QImg)defcreate_pixmap_tiles(self):"""
根据显示的图层创建pixmap的瓦片
:return:
"""
length =len(self.image.get_show_layer())
cur =0for key, value in self.image.get_show_layer().items():# 对于每个需要生成的瓦片单独创建一个pixmap
pixmap = self.create_pixmap(value)
self.pixmap_tiles[key]= pixmap
cur = cur +1
System.signal.signal_progress.emit("正在创建瓦片:", cur *100/ length)
QApplication.processEvents()print("pixmap瓦片创建完成")# self.image.start_image_static()
参考
[1] Python+GDAL栅格数据基本操作 https://blog.csdn.net/weixin_40625478/article/details/107839548
[2] Python空间数据处理1:GDAL读写遥感图像 https://blog.csdn.net/vonuo/article/details/74783291
[3]如何利用Python对遥感影像进行显示 https://blog.csdn.net/u012569919/article/details/109864492
[4] pyqt5 分割瓦片显示超大图像 https://blog.csdn.net/ZH_CS/article/details/124708547
版权归原作者 zstar-_ 所有, 如有侵权,请联系我们删除。