Linux下使用GDAL处理空间数据-猿码集

1. GDAL简介

GDAL (Geospatial Data Abstraction Library)是一组用于处理和转换各种空间数据格式的开源库。它提供了处理栅格数据、矢量数据和影像数据的功能，可以用于读取、写入和转换不同格式的空间数据。在Linux系统上，使用GDAL可以进行各种空间数据处理任务，包括空间数据格式转换、空间数据投影转换、空间数据裁剪、空间数据统计等。

2. 安装GDAL

在Linux系统上安装GDAL非常简单。可以通过系统的包管理器直接安装GDAL。

2.1 Ubuntu系统

在Ubuntu系统上，可以使用以下命令安装GDAL：

sudo apt-get install gdal-bin

2.2 CentOS系统

在CentOS系统上，可以使用以下命令安装GDAL：

sudo yum install gdal

2.3 编译安装

如果系统包管理器中没有GDAL的安装包，或者需要安装最新版本的GDAL，可以选择通过源码编译安装。

首先，从GDAL官方网站（https://gdal.org/）下载最新的GDAL源码压缩包。然后解压压缩包，并进入解压后的目录：

tar -zxvf gdal-x.x.x.tar.gz

cd gdal-x.x.x

在进入GDAL源码目录后，可以使用以下命令进行编译和安装：

./configure make

sudo make install

3. 使用GDAL处理空间数据

3.1 空间数据格式转换

GDAL可以将不同格式的空间数据相互转换，使得不同格式的空间数据可以在不同的GIS软件中进行使用。以下是一个示例代码，将Shapefile格式的矢量数据转换为GeoJSON格式：

import gdal
# 输入文件路径
input_file = "input.shp"
# 输出文件路径
output_file = "output.geojson"
# 打开输入文件
input_dataset = gdal.OpenEx(input_file)
# 创建输出驱动
output_driver = gdal.GetDriverByName("GeoJSON")
# 创建输出文件
output_dataset = output_driver.Create(output_file, 0, 0, 0, gdal.GDT_Unknown)
# 将输入文件中的图层写入输出文件
output_dataset.CopyLayer(input_dataset.GetLayer(), "")
# 关闭文件
input_dataset = None
output_dataset = None

在上述代码中，我们首先通过gdal.OpenEx函数打开输入文件，然后通过gdal.GetDriverByName函数创建输出驱动，接着使用output_driver.Create函数创建输出文件，并使用output_dataset.CopyLayer函数将输入文件中的图层写入输出文件。最后，通过将输入文件和输出文件设置为None来关闭文件。

3.2 空间数据投影转换

GDAL可以将空间数据的投影进行转换，使得不同投影的数据可以在同一个地理空间上进行叠加和分析。以下是一个示例代码，将一个Shapefile文件从WGS84投影转换为Mercator投影：

import gdal
# 输入文件路径
input_file = "input.shp"
# 输出文件路径
output_file = "output.shp"
# 打开输入文件
input_dataset = gdal.OpenEx(input_file)
# 创建输出驱动
output_driver = gdal.GetDriverByName("ESRI Shapefile")
# 创建输出文件
output_dataset = output_driver.Create(output_file, 0, 0, 0, gdal.GDT_Unknown)
# 获取输入文件的坐标转换对象
input_srs = input_dataset.GetLayer().GetSpatialRef()
# 创建输出文件的坐标转换对象
output_srs = gdal.osr.SpatialReference()
output_srs.ImportFromEPSG(3857)
# 创建坐标转换对象
transform = gdal.osr.CoordinateTransformation(input_srs, output_srs)
# 将输入文件中的图层投影转换后写入输出文件
gdal.VectorTranslate(output_dataset, input_dataset, options = ["-s_srs", input_srs.ExportToProj4(), "-t_srs", output_srs.ExportToProj4(), "-f", "ESRI Shapefile"])
# 关闭文件
input_dataset = None
output_dataset = None

在上述代码中，我们首先通过gdal.OpenEx函数打开输入文件，然后通过gdal.GetDriverByName函数创建输出驱动，接着使用output_driver.Create函数创建输出文件。通过input_dataset.GetLayer().GetSpatialRef()可以获取输入文件的坐标转换对象，通过gdal.osr.SpatialReference函数创建输出文件的坐标转换对象，然后使用gdal.osr.CoordinateTransformation函数创建坐标转换对象。最后，通过gdal.VectorTranslate函数将输入文件中的图层投影转换后写入输出文件，并将输入文件和输出文件设置为None来关闭文件。

3.3 空间数据裁剪

GDAL可以将空间数据根据给定的范围进行裁剪，使得只保留指定范围内的数据。以下是一个示例代码，裁剪一个GeoJSON文件，只保留指定范围内的数据：

import gdal
# 输入文件路径
input_file = "input.geojson"
# 输出文件路径
output_file = "output.geojson"
# 裁剪范围
xmin, ymin, xmax, ymax = 0, 0, 10, 10
# 打开输入文件
input_dataset = gdal.OpenEx(input_file)
# 创建输出驱动
output_driver = gdal.GetDriverByName("GeoJSON")
# 创建输出文件
output_dataset = output_driver.Create(output_file, 0, 0, 0, gdal.GDT_Unknown)
# 根据裁剪范围创建裁剪区域
clip_geom = gdal.osr.Geometry(gdal.osr.wkbLinearRing)
clip_geom.AddPoint(xmin, ymin)
clip_geom.AddPoint(xmin, ymax)
clip_geom.AddPoint(xmax, ymax)
clip_geom.AddPoint(xmax, ymin)
clip_geom.AddPoint(xmin, ymin)
clip_geom.CloseRings()
clip_poly = gdal.osr.Geometry(gdal.osr.wkbPolygon)
clip_poly.AddGeometry(clip_geom)
# 设置裁剪选项
options = gdal.VectorTranslateOptions(format="GeoJSON", clipSrc = clip_poly.ExportToWkt(), dstSRS='EPSG:4326')
# 执行裁剪操作
gdal.VectorTranslate(output_dataset, input_dataset, options=options)
# 关闭文件
input_dataset = None
output_dataset = None

在上述代码中，我们首先通过gdal.OpenEx函数打开输入文件，然后通过gdal.GetDriverByName函数创建输出驱动，接着使用output_driver.Create函数创建输出文件。通过gdal.osr.Geometry和gdal.osr.wkbLinearRing等类和函数可以创建裁剪区域的几何对象，然后使用gdal.osr.Geometry的ExportToWkt方法将裁剪区域导出为WKT格式。通过gdal.VectorTranslateOptions函数可以设置裁剪选项，然后通过gdal.VectorTranslate函数执行裁剪操作。

3.4 空间数据统计

GDAL可以计算空间数据的统计信息，包括最小值、最大值、平均值和标准差等。以下是一个示例代码，计算一个栅格数据的最小值、最大值和平均值：

import gdal
# 输入文件路径
input_file = "input.tif"
# 打开输入文件
input_dataset = gdal.OpenEx(input_file)
# 获取栅格数据的统计信息
min_value = input_dataset.GetRasterBand(1).GetMinimum()
max_value = input_dataset.GetRasterBand(1).GetMaximum()
mean_value = input_dataset.GetRasterBand(1).GetMean()
# 关闭文件
input_dataset = None
# 打印统计信息
print("Minimum value:", min_value)
print("Maximum value:", max_value)
print("Mean value:", mean_value)

在上述代码中，我们首先通过gdal.OpenEx函数打开输入文件，然后通过input_dataset.GetRasterBand(1)获取栅格数据的第一个波段，并使用GetMinimum、GetMaximum和GetMean方法获取统计信息。最后，通过将输入文件设置为None来关闭文件，并将统计信息打印出来。

总结

本文介绍了如何在Linux下使用GDAL处理空间数据。首先介绍了GDAL的简介和安装方法，然后详细介绍了GDAL的空间数据处理功能，包括空间数据格式转换、空间数据投影转换、空间数据裁剪和空间数据统计。通过示例代码和详细说明，帮助读者了解如何使用GDAL处理空间数据。

Linux下使用GDAL处理空间数据