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

这篇具有很好参考价值的文章主要介绍了【Python】GDAL基本操作/遥感大图显示。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

前言

遥感图像往往尺寸较大,无法用默认的图像浏览器加载。
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]采用了该方案进行图像显示。经实测,该方案能够有效解决遥感大图显示问题,并且拖动浏览较为流畅,但在显示之前需要耗费一定时间来切分瓦片。下面是瓦片显示实现的核心代码。

  1. 切分瓦片
    第一步是设定瓦片尺寸为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()
  1. 影像分析
    影像分析这步主要是用来统计更新中的影像极值,以便后续拉伸处理。
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 < 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()
  1. 初始化影像显示
    此步主要是对影像进行拉伸变换,将影像值拉伸到0-255之间,然后将处理后的数据保存进一个新的字典tiles_dict_show
def stretch_extremum(self, array):
    """
    影像拉伸-将影像值拉伸到0-255之间
    :param array:
    :return:
    """
    show_tiles = np.empty(array.shape, dtype='uint8')
    for band in range(self.bands):
        # 影像极值
        max_value = self.extremum_array[band][0]
        min_value = self.extremum_array[band][1]
        coefficient = (max_value - min_value) / 255
        if 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)
       
def init_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 = 0
    for 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()
  1. 创建瓦片并显示
    最后这步是根据每一个瓦片的数据,单独创建一个pixmap,从而能够使其在QGraphicsView进行显示。
def create_pixmap(self, image_array):
    """
    根据给定的数组生成pixmap
    :param image_array:
    :return:
    """
    shape = image_array.shape
    height = 0
    width = 0
    bytesPerComponent = 1
    if len(shape) == 2:
        height, width = shape
    elif len(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:
        # 如果不符合就不显示
        return
    return QPixmap.fromImage(QImg)
    
def create_pixmap_tiles(self):
    """
    根据显示的图层创建pixmap的瓦片
    :return:
    """
    length = len(self.image.get_show_layer())
    cur = 0
    for 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文章来源地址https://www.toymoban.com/news/detail-786148.html

到了这里,关于【Python】GDAL基本操作/遥感大图显示的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • Python 基本操作

    Python语法可以直接在命令行中编写并执行: 或者,可以在服务器上创建一个带有.py文件扩展名的Python文件,并在命令行中运行它: 缩进是指代码行开头的空格。 在其他编程语言中,代码中的缩进仅用于提高可读性,而在Python中,缩进非常重要。 Python使用缩进来表示一块代码

    2024年02月08日
    浏览(43)
  • Python基本操作

    啦啦啦,现在开始,打算做一期Python基础教程,欢迎大家来看哦! 这期文章真的是Python基础中的基础,相信有一定编程基础的小伙伴们都一定能看懂的… 本文共分为以下几个部分: 数与运算符 基本输入输出 注释 模块基本操作 小彩蛋 本文内容导图 — 数与运算符 算数运算符

    2024年02月11日
    浏览(38)
  • Python文件的基本操作

    文件的基本操作 文件的读写模式 文件的读写操作相关的方法 文件的操作模式 文件的练习题 文件的操作模式

    2024年02月10日
    浏览(37)
  • Opencv+Python图像基本操作

    目录 图像的读取、显示和保存 获取图像属性  图像截取  绘图功能 画线 画矩形 画圆圈 画椭圆          画多边形 向图像添加文本 cv2.imread() ,  cv2.imshow() ,  cv2.imwrite()分别表示读取图片,显示图片,写入图片   retval = cv2.imread(文件名 [,显示控制参数]) cv2.IMREAD_UNCHANGED:不改

    2024年02月04日
    浏览(51)
  • Python调用Gurobi基本操作

    接上篇学会了如何用python调用gurobipy之后,这篇总结一些学到的基本操作。 tuplelist、tupledict、multidict、创建list、 tuplelist是Python list的扩展对象,使用tuplelist()不能忘记from gurobipy import *,tuplelist增加了快速筛选select功能,比传统的if...else...筛选速度快。 tuplelist运行结果: for

    2023年04月25日
    浏览(36)
  • Python遥感图像处理应用篇036:GDAL+Scikit-image计算遥感图像梯度直方图HOG

    方向梯度直方图HOG(Histogram of Oriented Gradient, HOG)特征是一种在计算机视觉和图像处理中用来进行物体检测的特征描述子。它通过计算和统计图像局部区域的梯度方向直方图来构成特征。Hog特征结合SVM分类器已经被广泛应用于图像识别中,尤其在行人检测中获得了极大的成功

    2024年02月11日
    浏览(47)
  • Python中对基本文件操作

    保存数据放在磁盘中 f=open(‘文件’,‘w’)或者f=open(‘文件’,‘r’) 3.1 写数据(write) 如果文件不存在那么创建,如果存在那么就先清空,然后写入数据 对象=open(“文件”,w) 对象.write(“写入数据”) 对象.close 3.2读数据(read) 以读的方式打开文件,读取数据 对象 = open(“文

    2024年02月12日
    浏览(44)
  • Python中元组的基本操作

    元组(tuple)与列表类似,也是由一系列按特定顺序排列的元素组成,但与列表不同的是它是不可变序列。 定义元组时,可以创建一个空元组,也可以创建有值元组;如果元组中只有一个元素,元素后仍需要加英文逗号。另外元组中的元素可以是不同的数据类型。 输出以上

    2024年02月16日
    浏览(37)
  • Python 进阶(七): Word 基本操作

    Word 是一个十分常用的文字处理工具,通常我们都是手动来操作它,本节我们来看一下如何通过 Python 来操作。 Python 提供了   python-docx   库,该库就是为 Word 文档量身定制的,安装使用   pip install python-docx   命令即可。 首先,我们使用 Python 来创建一个 Word 文档并向其中写

    2024年02月08日
    浏览(44)
  • Python四种基本结构的操作

    SampleList = [] SampleList = list() append(obj) :在列表末尾添加元素obj extend(seq) :在列表末尾添加多个值,使用extend()函数,seq是一个可迭代对象,否则报错。 Insert(index,obj) :将元素obj添加到列表的index位置处。 Sample = [1,2,3,4,5] Sample.append(1) #[1,2,3,4,5,6] s = [7,8] Sample.extend(s) #[1,2,3,4,5,

    2024年02月06日
    浏览(37)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包