PythonGIS,GDAL# 读取栅格数据

PythonGIS,GDAL# 读取栅格数据

前言

新手上路，把我在学习pythongis中的过程记录，包括出现的问题以及解决办法记录下来，从而帮助更多的新手。我比较喜欢在Geany上操作，以下代码 也是在Geany上操作的。

NO.1 导入gdal、numpy包

（关于包的下载网址和装包过程大家可以搜索，网上很多教程）

from osgeo import gdal
#导入gdal包，这步骤也可验证你gdal包是否装好，环境设置是否正确。
import numpy as np
#导入numpy包，读取栅格数据的数值，数组矩阵运算

NO.2 读取栅格数据信息

#打开文件
#1 在cmd中可以cd到数据的文件夹，然后直接读取，cd的方法
[https://jingyan.baidu.com/article/7c6fb4282cd76880642c90b9.html]

dataset=gdal.Open(“XXX.tif”)
#2 在Geany中我是先把数据路径赋给fn，如下操作：
fn = r'D:\study\pythongis\chapter_02\XXX.tif'
dataset=gdal.Open(fn)#1 2方法选1种即可

#读取栅格矩阵的列、行、波段数（按顺序对应代码）
im_width=dataset.RasterXSize
im_height=dataset.RasterYSize
im_bands=dataset.RasterCount

#仿射矩阵，左上角像素的大地坐标和像素分辨率
im_geotrans=dataset.GetGeoTransform()
#地图投影信息，字符串表示
im_proj=dataset.GetProjection()
###这里我就出现了一些问题，error1-n:can't find proj.
###我在确保自己的gdal包安装正确且无丢失情况下，通过设置用户变量环境解决了。
###link:[关于GDAL运行时出现Cannot find proj.db的解决办法]
(https://blog.csdn.net/qq_31793023/article/details/103622134)
###大家可以尝试，若不懂该方法可以和我探讨。

#读取波段，参数为波段的索引号，波段索引号从1开始
band=dataset.GetRasterBand(1)
#用ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)，读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。
im_datas=band.ReadAsArray(0,0,im_width,im_height)

NO.3 读取栅格数据的值

#获取某个范围内（行列）的像元值；如1~5行和5~11列；（滤波计算中常用的）
data=im_datas[1:6,5:12]
print(data)
#输出结果：
      [[927, 934, 936, 938, 940, 946],
       [935, 950, 955, 959, 961, 965],
       [936, 950, 955, 959, 961, 965],
       [936, 959, 959, 961, 961, 966],
       [931, 958, 965, 962, 957, 950]]
del dataset#释放dataset,否则dataset会一直被占用，在arcgis或envi中打开该图像时显示文件已被占用

完成！
接下来会分享写入栅格数据的经验，然后是矢量数据的读写！