「GIS教程」基于Python

您所在的位置：网站首页 › gdal读取tif压缩类型 › 「GIS教程」基于Python

「GIS教程」基于Python

2023-03-21 17:36| 来源: 网络整理| 查看: 265

「GIS教程」基于Python-GDAL的实践总结发布时间: 2023-03-14 所属分类: GIS探秘文章标签: ArcPyGDAL 当前位置：麻辣GIS » GIS探秘 » 「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处理栅格数据，常用的属性与方法如下图所示。