「GIS教程」基于Python-GDAL的实践总结

1 介绍

最近对Python的gdal库比较感兴趣，汇总了一些常用代码，在这里进行下总结与分享。GDAL(Geospatial Data Abstraction Library) 是处理地理栅格数据和矢量数据的基础库，是开源地理空间基金会(Open Source Geospatial Foundation，OSGeo)的一个项目，其底层由C/C++实现，封装了Python、Java等语言的调用接口。GDAL具备多种栅格和矢量数据的读写、转换、处理能力，在Python中调用GDAL的API函数时，底层执行的是C/C++编译的二进制文件。

GDAL官网地址 https://gdal.org/

GDAL官网python-api文档 https://gdal.org/api/index.html#python-api

2 开发准备

这里为了方便演示测试，我们利用conda新建一个gdal的虚拟环境(用 Miniconda管理Python虚拟环境，简直不要太爽~)，安装大约需要三分钟。

# 创建一个名为gdal-test的python版本为3.8的虚拟环境
conda create -n gdal-test python=3.8 -c https://mirrors.aliyun.com/anaconda/pkgs/main/ -y
# 激活gdal-test环境
conda activate gdal-test
# 安装gdal,版本3.0.2
conda install gdal==3.0.2  -c https://mirrors.aliyun.com/anaconda/pkgs/main/ -y
# 安装matplotlib,版本3.5.0
conda install matplotlib==3.5.0  -c https://mirrors.aliyun.com/anaconda/pkgs/main/ -y
# 如果环境捣鼓坏了，可以删除了重新创建，删除命令是conda env remove -n gdal-test

如果全部执行完，如果运行程序提示找不到DLLL的话， 可从下面网站上下载GDAL的wheel文件进行安装。

加州大学的python扩展库 https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal

# 卸载gdal库
conda uninstall gdal
# 安装wheel文件,需要定位到文件的下载位置(例如：我下载到了本地G:\temp目录下)
pip install G:\temp\GDAL‑3.3.3‑cp38‑cp38‑win_amd64.whl

3 入门程序

3.1 读取TIF影像

• 效果预览
  

• 代码实现

# -*- coding: utf-8 -*-

# 导入gdal,注意导入的名称
from matplotlib import pyplot as plt # 画图模块
from osgeo import gdal  # 或者直接用import gdal

raster_path = r'./data-use/tif/AP_05726_FBS_F0680_RT1.dem.tif'

# 读取影像
dataset = gdal.Open(raster_path)
# 打开波段1（注意:用索引1，而不是0，来获取第一个波段）
band = dataset.GetRasterBand(1)
elevation = band.ReadAsArray()
print(elevation.shape)  # 行列数
print(elevation)  # 行列矩阵
# cmap='gist_earth'以gist_earth色带形式展示
plt.imshow(elevation, cmap='gist_earth')
plt.show()

3.2 读取HDF影像

• 效果预览
  

• 代码实现

# -*- coding: utf-8 -*-

import matplotlib.pylab as plt  # 画图模块
from osgeo import gdal, gdalconst

raster_path = r'G:\temp\MOSAIC_TMP_2019249.hdf'  # 输入你的hdf文件(绝对路径)

dataset = gdal.Open(raster_path, gdalconst.GA_ReadOnly)  # 读取遥感数据信息
band = dataset.GetRasterBand(1)
location_array = band.ReadAsArray()  # 二维行列像元矩阵
plt.title('gdal-hdf-show')  # 设置图框标题为gdal-hdf-show
plt.imshow(location_array, cmap='rainbow')  # 绘图，配色为rainbow
plt.colorbar()  # 绘制图例
plt.savefig("./results/20220917.04.png", dpi=300)
plt.show()

3.3 读取影像2.5维展示

• 效果预览
  

• 代码实现

# -*- coding: utf-8 -*-

import matplotlib.pylab as plt  # 画图模块
import numpy as np
from matplotlib import cm
from matplotlib.colors import LightSource
from osgeo import gdal

raster_path = r'./data-use/tif/AP_05726_FBS_F0680_RT1.dem.tif'
dataset = gdal.Open(raster_path)
adfGeoTransform = dataset.GetGeoTransform()  # 获取投影信息
band = dataset.GetRasterBand(1)  # 用gdal去读写你的数据，当然dem只有一个波段

nrows = dataset.RasterXSize
ncols = dataset.RasterYSize  # 这两个行就是读取数据的行列数

# 数据的平面四至
lonmin = adfGeoTransform[0]
latmin = adfGeoTransform[3]
lonmax = adfGeoTransform[0] + nrows * adfGeoTransform[1] + ncols * adfGeoTransform[2]
latmax = adfGeoTransform[3] + nrows * adfGeoTransform[4] + ncols * adfGeoTransform[5]

x = np.linspace(lonmin, lonmax, ncols)
y = np.linspace(latmin, latmax, nrows)

# 将数据的x，y，z化作numpy矩阵
lon, lat = np.meshgrid(x, y)
elevation = band.ReadAsArray(0, 0, nrows, ncols)
# 限定一个范围
region = np.s_[10:400, 10:400]
lon, lat, elevation = lon[region], lat[region], elevation[region]

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(12, 10))
ls = LightSource(270, 60)  # 设置你可视化数据的色带
rgb = ls.shade(elevation,
               cmap=cm.gist_earth,  # 设置颜色映射
               vert_exag=0.1,
               blend_mode='soft')
surf = ax.plot_surface(lon, lat, elevation,
                       rstride=1,  # 指定行的跨度
                       cstride=1,  # 指定列的跨度
                       facecolors=rgb,
                       linewidth=0,
                       antialiased=False, shade=False)
# 设置标题
plt.title("Python-gdal DEM show ", fontsize='large', fontweight='bold', color='#6666FF')
# fig.colorbar(surf, shrink=0.5, aspect=5)  # shrink越小，表示colorbar越小
plt.show()  # 最后渲染出2.5维图

3.4 读取栅格元数据

打开一幅栅格影像后，我们通过数据集+点的形式发现主要有如下六个方法或属性。比如说常见的投影信息，平面四至，像元数等等都可以获取到。

# -*- coding: utf-8 -*-
import os
from osgeo import gdal, osr

def get_geo_info(FileName):
    if os.path.exists(FileName) is False:
        raise Exception('[Errno 2] 该文件不存在: \'' + FileName + '\'')
    dataset = gdal.Open(FileName, gdal.GA_ReadOnly)
    if dataset == None:
        raise Exception('[Errno 3] 无法读取该文件: \'' + FileName + '\'')

    no_data_value = dataset.GetRasterBand(1).GetNoDataValue()  # 获取无效像元值
    xsize = dataset.RasterXSize  # 读取图像的宽度,x方向上的像素个数
    ysize = dataset.RasterYSize  # 读取图像的高度,y方向上的像素个数
    raster_count = int(dataset.RasterCount)  # 总像元数
    print('数据的大小（行/height，列/weight, RasterCount）：')
    print('(%s %s %s)' % (xsize, ysize, raster_count))
    geo_transform = dataset.GetGeoTransform()  # 获取
    projection = osr.SpatialReference()
    projection.ImportFromWkt(dataset.GetProjectionRef())
    data_type = dataset.GetRasterBand(1).DataType
    data_type = gdal.GetDataTypeName(data_type)
    # print('数据的类型{}'.format(data_type))
    return no_data_value, xsize, ysize, geo_transform, projection, data_type


if __name__ == '__main__':
    raster_path = r'./data-use/tif/NDVI_201405_bm.tif'
    info = get_geo_info(raster_path)
    for k in info:
        print(k)

通过一段时间的了解，发现gdal处理栅格数据，常用的属性与方法如下图所示。

4 问题解答

1. 为什么会写Python-GDAL这个专题？

   作为一个giser，总觉得打开桌面GIS软件，通过各种各样的工具，输入文件与参数执行，忽略了很多细节，这些细节也是GIS的魅力之一，通过编写脚本使我也更加理解了工具为什么要这么设计；再加上工作中有编写数据批处理脚本的需要，所以就总结了一些常用的脚本。

2. 为什么看了上面的内容后，感觉丝毫没有用处？

   法拉第曾回答:“新生的婴儿有什么用?新生的婴儿是会长大的”。通过练习与积累，遇到具体的问题，解决问题，这就是能力的提升。通过该专题的学习，你将更熟悉Python的基础语法，领悟它的魅力、巩固GIS基础知识，如OGC、格式转换、投影变化等等...最终达到通过模仿进行批处理脚本编写。

3. 为什么不使用ArcPy库，而利用GDAL库？

   arcpy也是很棒的一款空间数据处理库，工具全，教程多，操作简单方便，但它貌似不能直接pip安装，需要安装ArcGIS (pro)的产品才行，本着拥抱开源的态度，因此没有采用。

4. 后续会分享什么类型的脚本？

   目前打算先分享一些栅格数据处理的脚本，如数据转换（NC转TIF、ASC转TIF等）、投影转换、栅格处理（裁剪、重采样、NDVI计算、影像去云等等）

   

5 写在最后

以上就是本文的全部内容，代码和示例数据已上传至码云仓库，有需要自行下载查看，欢迎大家交流讨论~

https://gitee.com/fungiser/python-gdal-test