【Python】GDAL基本操作/遥感大图显示

#【Python】GDAL基本操作/遥感大图显示| 来源: 网络整理| 查看: 265

前言

遥感图像往往尺寸较大，无法用默认的图像浏览器加载。 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)来读取区域栅格图像数据，然后保存进字典。 class Tiles: 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 in range(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 in range(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 in range(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() 影像分析影像分析这步主要是用来统计更新中的影像极值，以便后续拉伸处理。 def image_extremum(self): length = len(self.tiles_dict_source) cur = 0 for key, value in self.tiles_dict_source.items(): # 统计极值 for band in range(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

【本文地址】

公司简介

联系我们