使用Python生成高程等值线图的DEM数据处理指南

本文还有配套的精品资源，点击获取

简介：在地理信息系统和遥感领域，数字高程模型（DEM）数据是关键，它们提供地形表面的连续高度信息，可用于多种应用。Python与GDAL/OGR、Rasterio和matplotlib等库结合，提供了一种高效处理DEM数据并生成高程等值线图的方法。该指南将指导用户理解DEM数据，掌握利用Python脚本处理这些数据，并通过可视化展示地形的起伏和高度变化。

1. 数字高程模型（DEM）介绍

数字高程模型（DEM）是一种在计算机中表示地形表面的方法，它通过一个规则的格网系统记录地面点的高程信息，被广泛应用于地形分析、地图制作和环境模拟等领域。DEM不仅能够提供精准的地形数据，还能通过各种方式分析地形特征，如坡度、坡向以及地形起伏度等。随着技术的发展，DEM数据也在不断地更新和优化，为地理信息系统（GIS）提供了强大的支持。在深入探讨如何使用Python处理DEM数据之前，了解DEM的基本概念及其应用场景是至关重要的。本章将从DEM的定义、数据结构以及在GIS中的应用等多个角度，为读者提供一个全面的DEM知识概览。

2. Python在GIS数据处理中的应用

2.1 Python在GIS数据处理中的优势与作用

Python已经成为数据科学、机器学习、网络开发等多个领域的首选语言，其在GIS数据处理中的应用也越来越广泛。原因在于Python具备一系列独特的语言特性和强大的生态系统。

2.1.1 Python语言特性及其在GIS中的应用前景

Python的简洁语法，使得代码易于编写和维护，非常适合快速开发。而其丰富的库支持，特别是GIS领域的专用库如GDAL/OGR、Shapely和Fiona等，为处理各种复杂的GIS数据提供了可能。Python的交互式特性使其非常适合进行数据分析和原型设计。

# 示例：使用Python和Shapely库来创建和操作几何对象
from shapely.geometry import Point, Polygon

point = Point(1, 2)
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

print(point.within(polygon))  # 检查点是否在多边形内

代码逻辑解读： 1. 导入Shapely库中的Point和Polygon类。 2. 创建一个点（1,2）和一个矩形多边形。 3. 使用 within 方法来检查点是否位于多边形内。 4. 执行后，输出布尔值表示点是否在多边形内。

Python的这些特性让它成为进行GIS数据处理的首选工具，尤其是在原型开发、数据清洗和复杂的空间分析等方面。

2.1.2 开源Python库在GIS数据处理中的角色

开源Python库为GIS数据处理提供了巨大的便利。GDAL/OGR是处理栅格和矢量数据的标准库，提供了读写大多数GIS数据格式的功能。此外，Pandas库可以用来处理表格数据，而NumPy库则在数值计算中扮演重要角色。通过这些库的组合使用，可以实现从数据读取到分析再到可视化的完整GIS数据处理流程。

import pandas as pd
import numpy as np

# 示例：使用Pandas读取CSV文件中的GIS表格数据
data = pd.read_csv('gis_data.csv')

# 数据处理
data = data[['longitude', 'latitude']]
data = data.dropna()

# 使用NumPy进行计算
mean_lat = np.mean(data['latitude'])
print(f'平均纬度是 {mean_lat}')

代码逻辑解读： 1. 导入Pandas和NumPy库。 2. 读取存储GIS数据的CSV文件到DataFrame。 3. 删除数据中的空值。 4. 计算纬度的平均值并打印结果。

这些开源库不仅简化了数据处理流程，而且在很大程度上提高了开发效率和代码的可读性。

2.2 Python进行DEM数据预处理

数字高程模型（DEM）是一种重要的地理信息数据，可以表示地表的高程信息。Python在处理DEM数据方面具有强大的数据处理能力。

2.2.1 读取DEM数据与初步分析

在进行DEM数据预处理之前，需要先读取数据。Python中GDAL库的gdal.Open方法是一个常用的函数，可以用来打开栅格数据文件。

from osgeo import gdal

# 使用GDAL打开DEM数据文件
dataset = gdal.Open('dem_file.tif')

# 获取栅格数据集的地理转换信息
geo_transform = dataset.GetGeoTransform()
print(geo_transform)

代码逻辑解读： 1. 导入GDAL模块。 2. 使用gdal.Open函数打开名为'dem_file.tif'的DEM数据文件。 3. 通过GetGeoTransform方法获取数据的地理变换参数。 4. 打印地理变换参数。

这个过程包括了数据读取和获取基础信息的初步分析。了解这些信息对于后续的数据处理至关重要。

2.2.2 数据清洗与格式转换技巧

DEM数据在使用之前往往需要进行数据清洗和格式转换。例如，有时需要将DEM数据转换为其他格式，以便与其他系统兼容。Python中的GDAL库提供了强大的格式转换能力，可以将DEM数据轻松转换为所需格式。

# 使用GDAL库将DEM数据转换为其他格式
driver = gdal.GetDriverByName('GTiff')
out_dataset = driver.CreateCopy('output_dem.tif', dataset, 0)
out_dataset = None

代码逻辑解读： 1. 使用GetDriverByName方法获取TIFF格式驱动。 2. 创建一个输出数据集，将原数据集复制到新文件。 3. 输出结果，关闭数据集。

通过这个简单的代码块，Python和GDAL完成了DEM数据从读取到输出的完整转换过程，展示了其在数据处理中的强大能力。

3. DEM数据存储格式（TIFF）

3.1 TIFF格式概述及其在DEM中的应用

3.1.1 TIFF格式的特点及优势分析

TIFF（Tagged Image File Format）是一种灵活的位图图像格式，它支持多种数据类型和图像压缩方式。在GIS领域，尤其是在数字高程模型（DEM）数据处理中，TIFF格式因其开放性、可扩展性和高效性被广泛应用。

灵活性

TIFF支持多种像素格式，包括灰度、RGB和CMYK等，以及不同位深度，如8位、16位、32位等，这使得TIFF非常适合存储高精度的DEM数据。

开放性

TIFF格式是开放的，不依赖于特定的软件或硬件平台。这意味着TIFF文件可以被多种GIS软件读取和处理，包括开源工具如QGIS和商业软件如ArcGIS。

可扩展性

TIFF文件可以包含多种“标签”，这些标签可以存储额外的信息，如图像的地理参考数据。这对于GIS应用而言极为重要，因为它允许将空间坐标直接嵌入到文件中，从而便于在DEM分析中进行空间定位。

高效性

TIFF格式支持无损压缩，这有助于减小文件大小，同时保持图像质量。对于存储大型DEM数据集来说，这意味着可以节省存储空间和传输时间，同时维持数据的完整性。

3.1.2 TIFF格式与DEM数据的关系

DEM数据通常由一系列规则的格网组成，每个格网单元对应一个高程值。TIFF格式能够很好地表达这种规则的栅格数据结构，使得DEM数据能够以一种标准化、高效的方式被存储和传输。

栅格数据表示

TIFF格式特别适合存储栅格数据，因此它可以无缝地表示DEM中的每个高程点。这使得TIFF成为存储大规模地形数据的理想选择。

地理参考信息

通过TIFF的标签系统，地理坐标和投影信息可以直接嵌入到DEM数据文件中。这样，DEM数据不仅包含了高程信息，还包含了定位信息，为后续的空间分析和地图制作提供了极大的便利。

多波段存储能力

TIFF支持多波段存储，这意味着一个TIFF文件可以同时存储多个图像通道，例如，在表示彩色高程图时，可以将高程信息分别存储在不同的颜色通道中。

3.2 TIFF格式数据的读写与管理

3.2.1 使用Python读写TIFF文件

Python作为一个功能强大的编程语言，在读写TIFF文件方面也表现得游刃有余。尤其是借助一些专门的库如 rasterio ，可以非常方便地处理TIFF文件。

安装必要的Python库

首先，需要安装 rasterio 库，可以使用pip进行安装：

pip install rasterio

基本读取操作

下面的Python代码展示了如何使用 rasterio 读取一个TIFF文件并打印其基本信息：

import rasterio

# 打开TIFF文件
with rasterio.open('dem_data.tif') as dataset:
    # 打印文件信息
    print(dataset.meta)

数据读取与分析

进一步，我们可以读取TIFF文件中的数据，并进行一些基础分析。例如，下面的代码读取了TIFF文件中的所有高程值，并计算平均高程：

import numpy as np

# 读取高程数据
with rasterio.open('dem_data.tif') as dataset:
    elevation_data = dataset.read(1)
    mean_elevation = np.mean(elevation_data)
    print(f'平均高程: {mean_elevation}')

3.2.2 TIFF文件的属性设置与数据检索

属性设置

TIFF文件的属性设置包括图像的分辨率、空间参考系统、波段信息等。通过 rasterio 库，我们可以设置或修改这些属性。例如，下面的代码展示了如何设置TIFF文件的地理变换参数：

with rasterio.open(
    'dem_data.tif', 
    'w',
    driver='GTiff',
    height=dataset.height,
    width=dataset.width,
    count=1,
    dtype=dataset.dtype,
    crs=dataset.crs,
    transform=dataset.transform,
) as dst:
    dst.write(elevation_data, 1)

数据检索

数据检索指的是根据地理坐标来检索TIFF文件中的具体像素值。 rasterio 库提供的方法可以轻松实现这一点。下面的代码段展示了如何根据坐标点检索其对应的高程值：

# 坐标点
x, y = 440000, 3751000
# 检索高程值
with rasterio.open('dem_data.tif') as dataset:
    row, col = dataset.index(x, y)
    elevation_value = dataset.read(1, window=rio.windows.Window(col, row, 1, 1))
    print(f'在坐标({x}, {y})处的高程值为: {elevation_value[0][0]}')

通过上述操作，我们不仅了解了如何读写TIFF文件，还学会了如何检索和分析TIFF格式的DEM数据。这些技能对于GIS数据处理来说至关重要，并可以进一步用于空间分析、地形建模等领域。

4. Python库GDAL/OGR的DEM数据处理功能

GDAL/OGR库是地理空间数据抽象库，广泛应用于遥感图像处理和地理信息系统(GIS)领域。GDAL主要用于栅格数据处理，而OGR专注于矢量数据。在DEM数据处理中，GDAL/OGR提供了强大的工具来读取、转换、分析和导出数据。

4.1 GDAL/OGR库的介绍与安装

4.1.1 GDAL/OGR库的基本概念与作用

GDAL (Geospatial Data Abstraction Library) 和 OGR (OpenGIS Simple Features for GIS) 是两个独立但相互关联的库。GDAL用于栅格数据处理，可以读取和写入多种栅格数据格式，而OGR专注于矢量数据格式的读写和转换。两者共同构成一个完整的地理空间数据抽象层，允许开发者不必关心数据格式的差异即可进行数据处理。

GDAL/OGR库的作用主要体现在以下几个方面：

1. 数据格式抽象 ：通过统一的API接口，简化了不同地理空间数据格式之间的转换工作。
2. 数据处理能力 ：提供了一系列命令行工具和编程接口，能够执行多种空间数据操作，如重投影、裁剪、像素值统计等。
3. 空间分析和处理 ：支持多种空间分析功能，如缓冲区分析、空间关系判断等。
4. 格式支持广泛 ：支持多种矢量和栅格格式，包括一些专有格式和一些较不常见的开源格式。

4.1.2 安装GDAL/OGR库的步骤与环境配置

安装GDAL/OGR库相对直接，以下是使用Python进行DEM数据处理的环境搭建步骤：

1. 安装Python环境 ：确保系统中已经安装了Python。GDAL/OGR提供了Python绑定，可以通过Python脚本使用GDAL/OGR的功能。

2. 安装GDAL/OGR库 ：可以通过包管理器安装GDAL/OGR，例如在Ubuntu系统中，可以使用以下命令：

bash sudo apt-get install python-gdal

或者使用Python的包管理工具pip：

bash pip install gdal

1. 环境配置 ：确保GDAL/OGR库的路径被添加到了系统的环境变量中，这样Python脚本才能正确加载库。如果使用的是Windows系统，需要在系统环境变量中添加GDAL/OGR的安装路径。

2. 验证安装 ：安装完成后，可以使用Python交互式环境或者创建一个Python脚本来导入GDAL/OGR库，以确认安装成功。

python from osgeo import gdal, ogr print(gdal.__version__) print(ogr.__version__)

如果上述代码可以正常执行并输出版本信息，说明GDAL/OGR库安装正确。

4.2 GDAL/OGR在DEM数据操作中的应用

4.2.1 DEM数据的读取与导出

GDAL库提供了读取和导出DEM数据的强大功能。DEM数据通常以栅格形式存储，利用GDAL可以轻松读取这些数据，并将它们转换为其他格式或者进行进一步的处理。

读取DEM数据

利用GDAL库读取DEM数据十分简单。以下是一个基本的Python脚本示例，展示如何读取一个DEM文件，并获取其基本信息：

from osgeo import gdal

# 打开DEM文件
dataset = gdal.Open('path_to_your_dem.tif')

# 获取栅格数据集的基本信息
print(f"Driver: {dataset.GetDriver().ShortName}/{dataset.GetDriver().LongName}")
print(f"Size is {dataset.RasterXSize}x{dataset.RasterYSize}x{dataset.RasterCount}")
print(f"Projection is `{dataset.GetProjection()}`")
geo_transform = dataset.GetGeoTransform()
if geo_transform:
    print(f"Origin = ({geo_transform[0]}, {geo_transform[3]})")
    print(f"Pixel Size = ({geo_transform[1]}, {geo_transform[5]})")

# 关闭数据集
dataset = None

导出DEM数据

在处理完DEM数据后，可能需要将数据导出为不同的格式。利用GDAL，可以轻松地将处理后的数据保存为其他格式：

from osgeo import gdal

# 打开原始DEM数据集
source_dataset = gdal.Open('path_to_your_dem.tif', gdal.GA_Update)

# 创建新的文件，并设置数据类型和大小
driver = gdal.GetDriverByName('GTiff')
target_dataset = driver.Create('path_to_your_output.tif', source_dataset.RasterXSize, source_dataset.RasterYSize, source_dataset.RasterCount, gdal.GDT_Float32)

# 复制数据集元数据
target_dataset.SetGeoTransform(source_dataset.GetGeoTransform())
target_dataset.SetProjection(source_dataset.GetProjection())

# 复制栅格数据
for i in range(1, source_dataset.RasterCount + 1):
    band = source_dataset.GetRasterBand(i)
    target_band = target_dataset.GetRasterBand(i)
    target_band.WriteArray(band.ReadAsArray())

# 清理并关闭数据集
target_dataset = None
source_dataset = None

以上代码段展示了如何将一个DEM数据集以TIFF格式保存到新文件中，同时也展示了如何设置数据类型和元数据。

4.2.2 DEM数据的空间分析与处理

GDAL/OGR不仅可以读取和导出数据，还提供了对DEM数据进行空间分析和处理的强大功能。这些功能对于地形分析、水文建模和可视化等应用尤为重要。

DEM数据重投影

重投影是空间数据处理中的常见需求。由于不同地理空间数据可能使用不同的坐标系，所以在进行空间分析或数据整合前需要将它们投影到统一的坐标系下。

以下是使用GDAL库进行DEM数据重投影的一个基本示例：

from osgeo import gdal

# 打开原始DEM数据集
dataset = gdal.Open('path_to_your_input_dem.tif')

# 获取源和目标坐标系
source_prj = dataset.GetProjection()
target_prj = "+proj=longlat +datum=WGS84"

# 获取仿射变换参数
geotransform = dataset.GetGeoTransform()

# 创建输出文件，并设置目标坐标系
driver = gdal.GetDriverByName('GTiff')
output_dataset = driver.Create('path_to_your_output_dem.tif', dataset.RasterXSize, dataset.RasterYSize, dataset.RasterCount, gdal.GDT_Float32)
output_dataset.SetProjection(target_prj)

# 设置输出文件的仿射变换参数
output_dataset.SetGeoTransform(geotransform)

# 使用重投影算法将数据从源数据集复制到输出数据集
for i in range(1, dataset.RasterCount + 1):
    input_band = dataset.GetRasterBand(i)
    output_band = output_dataset.GetRasterBand(i)
    output_band.SetNoDataValue(input_band.GetNoDataValue())
    output_band.WriteArray(input_band.ReadAsArray())

# 清理并关闭数据集
output_dataset = None
dataset = None

该代码示例展示了如何将DEM数据从一个坐标系投影转换到WGS84坐标系。

DEM数据裁剪

裁剪是另一种常见的DEM数据处理操作。在只需要分析特定区域的情况下，可以使用GDAL来裁剪DEM数据。

from osgeo import gdal

# 打开原始DEM数据集
dataset = gdal.Open('path_to_your_input_dem.tif')

# 设置裁剪区域
ul_x, ul_y = (100000, 200000)  # 裁剪区域左上角坐标
lr_x, lr_y = (200000, 100000)  # 裁剪区域右下角坐标

# 计算裁剪区域的大小
x_size = int(lr_x - ul_x)
y_size = int(lr_y - ul_y)

# 裁剪区域仿射变换参数
geotransform = (ul_x, dataset.GetGeoTransform()[1], 0, ul_y, 0, dataset.GetGeoTransform()[5])

# 创建输出文件
driver = gdal.GetDriverByName('GTiff')
clipped_dataset = driver.Create('path_to_your_clipped_dem.tif', x_size, y_size, 1, gdal.GDT_Float32)
clipped_dataset.SetGeoTransform(geotransform)
clipped_dataset.SetProjection(dataset.GetProjection())

# 从输入数据集中复制裁剪区域数据到输出数据集
band = dataset.GetRasterBand(1)
band.ReadAsArray(ul_x, ul_y, x_size, y_size, x_size, y_size, buf_type=gdal.GDT_Float32)
clipped_band = clipped_dataset.GetRasterBand(1)
clipped_band.WriteArray(band.ReadAsArray())

# 清理并关闭数据集
clipped_dataset = None
dataset = None

这个示例演示了如何对DEM数据进行裁剪操作。

GDAL/OGR库的DEM数据处理功能十分强大且具有灵活性，能够在GIS数据处理中提供极大的便利。通过掌握GDAL/OGR库，可以将DEM数据的读取、处理和分析变得更加高效和方便。

5. 使用Python绘制高程等值线图

在地理信息系统（GIS）的众多应用中，绘制高程等值线图是一项基础而重要的任务，这可以帮助我们直观地理解地形的起伏变化。Python作为一个强大的编程语言，配合其丰富的数据处理库，使得绘制高程等值线图变得简单而高效。

5.1 matplotlib库的基本使用方法

matplotlib是一个用于创建静态、动画和交互式可视化的Python库，它提供了绘图函数和对象层次结构，是进行科学计算可视化的重要工具。

5.1.1 matplotlib库在数据可视化中的作用

matplotlib可以绘制各种类型的图表，包括折线图、散点图、条形图、柱状图、直方图、饼图、等值线图等。它还能将图表保存为多种格式的图片文件，并可以与其它GUI工具包集成。在GIS应用中，matplotlib经常被用来将数据分析结果转换为可视化图表，帮助人们理解和解释数据。

5.1.2 基础绘图技能与等值线图制作技巧

在使用matplotlib绘图时，主要步骤包括准备数据、创建图表对象、绘制数据、添加图例和标签等。绘制等值线图时，通常需要X、Y坐标数据和对应的Z值，这些数据可以从DEM文件中获取。

下面是一个简单的等值线图绘制实例：

import numpy as np
import matplotlib.pyplot as plt

# 创建X和Y坐标网格
x = np.linspace(0, 10, 200)
y = np.linspace(0, 10, 200)
X, Y = np.meshgrid(x, y)

# 生成Z数据，这里为X+Y的函数值，实际中通常是DEM数据
Z = np.sin(np.sqrt(X**2 + Y**2))

# 绘制等值线图
plt.figure(figsize=(10, 6))
plt.contour(X, Y, Z, 20, cmap='RdGy')  # 20是等值线的数量，cmap是颜色映射

# 添加图例和标签
plt.colorbar(label='高程')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.title('高程等值线图')

# 显示图表
plt.show()

5.2 高程等值线图的应用实例

等值线图在地形分析、水文学和土木工程等领域中有着广泛应用。它们可以清晰地显示地形坡度、高度变化趋势和潜在的危险区域。

5.2.1 等值线图在地形分析中的重要性

通过等值线图，我们可以直观地看到地形的高程变化，比如山脊、山谷和盆地的位置，这有助于进行地形评估和理解地形对环境的影响。

5.2.2 利用matplotlib绘制高程等值线图的实战演练

下面我们通过实际的DEM数据来绘制高程等值线图。假设我们有一个TIFF格式的DEM文件，我们可以使用GDAL/OGR库来读取数据，并使用matplotlib来绘制等值线图。

from osgeo import gdal
import matplotlib.pyplot as plt

# 打开TIFF文件
dataset = gdal.Open('dem_data.tif')
band = dataset.GetRasterBand(1)

# 读取DEM数据
dem_array = band.ReadAsArray()

# 获取X和Y的起始坐标以及栅格的宽度和高度
geo_transform = dataset.GetGeoTransform()
x_min, x_size, _, y_max, _, y_size = geo_transform

# 生成X和Y坐标网格
x = np.linspace(0, x_size, dem_array.shape[1])
y = np.linspace(0, y_size, dem_array.shape[0])
X, Y = np.meshgrid(x, y)

# 绘制等值线图
plt.figure(figsize=(12, 8))
plt.contour(X, Y, dem_array, 20, cmap='viridis')

# 添加图例和标签
plt.colorbar(label='高程 (米)')
plt.xlabel('经度')
plt.ylabel('纬度')
plt.title('实际高程等值线图')

# 显示图表
plt.show()

# 关闭数据集
dataset = None

这段代码首先读取一个TIFF格式的DEM文件，然后使用 ReadAsArray 方法获取其中的高程数据。接着，利用读取的DEM数据和GeoTransform参数来创建X和Y坐标网格，并使用matplotlib的 contour 函数绘制等值线图。最后，添加必要的图例、标签，并展示图表。通过这样的流程，我们可以将DEM数据可视化为直观的等值线图，进而分析地形特征。

本文还有配套的精品资源，点击获取

简介：在地理信息系统和遥感领域，数字高程模型（DEM）数据是关键，它们提供地形表面的连续高度信息，可用于多种应用。Python与GDAL/OGR、Rasterio和matplotlib等库结合，提供了一种高效处理DEM数据并生成高程等值线图的方法。该指南将指导用户理解DEM数据，掌握利用Python脚本处理这些数据，并通过可视化展示地形的起伏和高度变化。

本文还有配套的精品资源，点击获取