Python引入GDAL进行DEM数据信息获取

Python进行DEM数据信息获取

实现思路

1. 采用gdal进行DEM（栅格）数据读取
2. 运用gdalJ进行DEM数据的信息获取

下载和引入GDAL包

可以有两种方式进行gdal包的安装：

• 直接用pip下载

  pip install gdal

  （该方法我安装不成功，显示需要安装Microsorft Visual C++ 14.0）

• 通过下载python非官方文件，下载gdal包

  1. python拓展包下载地址：[Unofficial Windows Binaries for Python Extension Packages](Python Extension Packages for Windows - Christoph Gohlke (uci.edu))

  2. 选择与自己python版本对应的whl包，其中cpXX-cpXX数字表示版本信息
     

  3. 打开命令行，路径转入到gdal.whl的路径下，运行：
     pip install GDAL‑3.3.1‑cpXX‑cpXX‑win_amd64.whl

引入gdal包：

from osgeo import gdal

DEM等栅格数据的信息获取

DEM概况信息（行列、坐标、投影信息等）

dataset = gdal.Open('./data/dem.tif') # 打开DEM数据

dem_XSize = dataset.RasterXSize  # 列数
dem_YSize = dataset.RasterYSize  # 行数
dem_bands = dataset.RasterCount  # 波段数

dem_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
dem_proj = dataset.GetProjection()  # 地图投影信息

band = dataset.GetRasterBand(1) # 获取第一个波段的信息

dem_data = band.ReadAsArray(0, 0, dem_XSize, dem_YSize) # 读取数据

读取每个单元的高程信息

for i in range(dem_YSize): #遍历行
    for j in range(dem_XSize): # 遍历列
        cur = dem_data[i][j] # dem矩阵的值（高程）

读取DEM坐标信息

def getLngLat(row, col, geotransform):
    px = geotransform[0] + col * geotransform[1] + row * geotransform[2]
    py = geotransform[3] + col * geotransform[4] + row * geotransform[5]
    return [px, py]
    
for i in range(dem_YSize): #遍历行
    for j in range(dem_XSize): # 遍历列
        pos = getLngLat(i,j,dem_geotrans) # 坐标

完整代码

# -*- coding: utf-8 -*-
import math
from osgeo import gdal
import numpy as np

dataset = gdal.Open('./data/dem.tif')

dem_XSize = dataset.RasterXSize  # 列数
dem_YSize = dataset.RasterYSize  # 行数
dem_bands = dataset.RasterCount  # 波段数
# print(dem_XSize, dem_YSize, dem_bands)

dem_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
# print("左上角地理坐标", dem_geotrans[0], dem_geotrans[3])
dem_proj = dataset.GetProjection()  # 地图投影信息

# 读取某一像素点的值
# 读取一个波段，其参数为波段的索引号，波段索引号从1开始
band = dataset.GetRasterBand(1)

# 用ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)，读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。以下为读取整幅图像
dem_data = band.ReadAsArray(0, 0, dem_XSize, dem_YSize)
# print(dem_data)

# 获取经纬度
def getLngLat(row, col, geotransform):
    px = geotransform[0] + col * geotransform[1] + row * geotransform[2]
    py = geotransform[3] + col * geotransform[4] + row * geotransform[5]
    return [px, py]

rowLen = dem_YSize
colLen = dem_XSize
rowStep = 1
colStep = 1
data = dem_data[0:rowLen, 0:colLen]
curRow = 0  # 当前行
while curRow < rowLen:
    curCol = 0  # 当前列
    while curCol < colLen:
        cur_pos = getLngLat(curRow, curCol, dem_geotrans) # 坐标
        cur_height = dem_data[curRow][curCol] # dem矩阵的值（高程）
        curCol += colStep
    curRow += rowStep
# 释放内存。如果不释放，在arcgis或envi中打开该图像时显示文件已被占用
del dataset

源码