Python库：gdal_IT分享知识网

大家好，欢迎来到IT知识分享网。

文章目录

GDAL（Geospatial Data Abstraction Library）是处理地理空间数据的开源库，广泛用于地理信息系统（GIS）和遥感领域。GDAL 提供了读取、写入、转换和处理栅格数据（如卫星图像、DEM等）和矢量数据（如Shapefile）的功能。

Python 通过 gdal 库提供了对 GDAL 的绑定，使得开发者可以在 Python 环境中使用 GDAL 的功能。

1. 安装 `gdal` 库

在安装 gdal 库之前，需要确保已经安装了 GDAL 库本身。然后可以通过 pip 安装 Python 的 gdal 绑定：

pip install gdal

2. 基本概念

在使用 gdal 库之前，了解一些基本概念是很有帮助的：

栅格数据：栅格数据是由像素组成的图像，每个像素都有一个值，通常表示某种测量值（如高度、温度、反射率等）。
矢量数据：矢量数据由点、线、面等几何对象组成，通常用于表示地理特征（如道路、建筑物、河流等）。
投影：地理数据通常需要投影到二维平面上，常见的投影方式包括 UTM、WGS84 等。
坐标系：地理数据通常使用特定的坐标系来表示位置，如 WGS84（经纬度）、UTM 等。

3. 读取栅格数据

gdal 库可以读取多种栅格数据格式，如 GeoTIFF、JPEG、PNG 等。以下是一个简单的示例，展示如何读取一个 GeoTIFF 文件并获取其基本信息。

from osgeo import gdal # 打开栅格数据文件 dataset = gdal.Open('example.tif') # 获取栅格数据的宽度和高度 width = dataset.RasterXSize height = dataset.RasterYSize # 获取栅格数据的波段数量 bands = dataset.RasterCount # 获取栅格数据的投影信息 projection = dataset.GetProjection() # 获取栅格数据的地理变换信息（仿射变换参数） geotransform = dataset.GetGeoTransform() # 获取第一个波段的数据 band = dataset.GetRasterBand(1) data = band.ReadAsArray() # 打印基本信息 print(f"Width: { 
      width}, Height: { 
      height}, Bands: { 
      bands}") print(f"Projection: { 
      projection}") print(f"GeoTransform: { 
      geotransform}") print(f"First band data shape: { 
      data.shape}")

4. 读取矢量数据

gdal 库也可以读取矢量数据，如 Shapefile。以下是一个简单的示例，展示如何读取一个 Shapefile 并获取其几何对象和属性。

from osgeo import ogr # 打开矢量数据文件 driver = ogr.GetDriverByName('ESRI Shapefile') dataSource = driver.Open('example.shp', 0) # 获取图层 layer = dataSource.GetLayer() # 遍历图层中的每个要素 for feature in layer: # 获取几何对象 geometry = feature.GetGeometryRef() print(f"Geometry: { 
      geometry.GetGeometryName()}") # 获取属性 attributes = feature.items() print(f"Attributes: { 
      attributes}")

5. 写入栅格数据

gdal 库不仅可以读取数据，还可以创建和写入栅格数据。以下是一个简单的示例，展示如何创建一个新的 GeoTIFF 文件并写入数据。

from osgeo import gdal, osr import numpy as np # 创建一个新的栅格数据集 driver = gdal.GetDriverByName('GTiff') dataset = driver.Create('output.tif', 512, 512, 1, gdal.GDT_Float32) # 设置地理变换参数 geotransform = (0, 1, 0, 0, 0, -1) dataset.SetGeoTransform(geotransform) # 设置投影信息 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # WGS84 dataset.SetProjection(srs.ExportToWkt()) # 创建一个波段并写入数据 band = dataset.GetRasterBand(1) data = np.random.rand(512, 512).astype(np.float32) band.WriteArray(data) # 关闭数据集 dataset = None

6. 写入矢量数据

同样，gdal 库也可以创建和写入矢量数据。以下是一个简单的示例，展示如何创建一个新的 Shapefile 并写入几何对象和属性。

from osgeo import ogr, osr # 创建一个新的矢量数据集 driver = ogr.GetDriverByName('ESRI Shapefile') dataSource = driver.CreateDataSource('output.shp') # 创建一个图层 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # WGS84 layer = dataSource.CreateLayer('output', srs, ogr.wkbPoint) # 创建一个字段 field_defn = ogr.FieldDefn('Name', ogr.OFTString) field_defn.SetWidth(24) layer.CreateField(field_defn) # 创建一个要素并写入几何对象和属性 feature_defn = layer.GetLayerDefn() feature = ogr.Feature(feature_defn) feature.SetField('Name', 'Point 1') # 创建几何对象 point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(10, 20) feature.SetGeometry(point) # 将要素写入图层 layer.CreateFeature(feature) # 关闭数据源 dataSource = None

7. 常见操作

7.1 重采样（Resampling）

重采样是指改变栅格数据的分辨率，即改变每个像素的大小。常见的重采样方法包括最近邻插值、双线性插值、三次卷积插值等。

from osgeo import gdal, gdalconst # 打开原始数据集 src_ds = gdal.Open('input.tif') # 设置目标分辨率 x_res = src_ds.RasterXSize * 2 # 例如，将宽度增加一倍 y_res = src_ds.RasterYSize * 2 # 例如，将高度增加一倍 # 创建目标数据集 dst_ds = gdal.GetDriverByName('GTiff').Create('output_resampled.tif', x_res, y_res, src_ds.RasterCount, gdal.GDT_Float32) # 设置地理变换参数 geotransform = list(src_ds.GetGeoTransform()) geotransform[1] /= 2 # 调整分辨率 geotransform[5] /= 2 # 调整分辨率 dst_ds.SetGeoTransform(geotransform) # 设置投影信息 dst_ds.SetProjection(src_ds.GetProjection()) # 执行重采样 gdal.ReprojectImage(src_ds, dst_ds, src_ds.GetProjection(), src_ds.GetProjection(), gdalconst.GRA_Bilinear) # 关闭数据集 src_ds = None dst_ds = None

7.2 裁剪（Clipping）

裁剪是指根据几何对象（如多边形）裁剪栅格数据，只保留几何对象内的部分。

from osgeo import gdal, ogr, osr # 打开原始数据集 src_ds = gdal.Open('input.tif') # 打开几何对象数据集（例如，一个Shapefile） shp_ds = ogr.Open('clip_polygon.shp') layer = shp_ds.GetLayer() # 创建目标数据集 dst_ds = gdal.GetDriverByName('GTiff').Create('output_clipped.tif', src_ds.RasterXSize, src_ds.RasterYSize, src_ds.RasterCount, gdal.GDT_Float32) # 设置地理变换参数和投影信息 dst_ds.SetGeoTransform(src_ds.GetGeoTransform()) dst_ds.SetProjection(src_ds.GetProjection()) # 执行裁剪 gdal.Warp(dst_ds, src_ds, cutlineDSName='clip_polygon.shp', cropToCutline=True) # 关闭数据集 src_ds = None dst_ds = None shp_ds = None

7.3 投影转换（Reprojection）

投影转换是指将数据从一个坐标系转换到另一个坐标系。常见的投影包括 WGS84、UTM 等。

from osgeo import gdal, osr # 打开原始数据集 src_ds = gdal.Open('input.tif') # 设置目标投影（例如，UTM Zone 18N） dst_srs = osr.SpatialReference() dst_srs.ImportFromEPSG(32618) # UTM Zone 18N # 创建目标数据集 dst_ds = gdal.GetDriverByName('GTiff').Create('output_reprojected.tif', src_ds.RasterXSize, src_ds.RasterYSize, src_ds.RasterCount, gdal.GDT_Float32) # 执行投影转换 gdal.Warp(dst_ds, src_ds, dstSRS=dst_srs.ExportToWkt()) # 关闭数据集 src_ds = None dst_ds = None

7.4 数据格式转换（Format Conversion）

数据格式转换是指将数据从一种格式转换为另一种格式。常见的栅格数据格式包括 GeoTIFF、JPEG、PNG 等。

from osgeo import gdal # 打开原始数据集 src_ds = gdal.Open('input.tif') # 创建目标数据集（例如，转换为JPEG格式） dst_ds = gdal.GetDriverByName('JPEG').CreateCopy('output_converted.jpg', src_ds) # 关闭数据集 src_ds = None dst_ds = None

8. 总结

gdal 库是一个功能强大的工具，适用于处理地理空间数据。通过 Python 的 gdal 绑定，开发者可以在 Python 环境中轻松地读取、写入、转换和处理栅格和矢量数据。无论是进行简单的数据读取，还是复杂的图像处理和分析，gdal 库都能提供强大的支持。

免责声明：本站所有文章内容,图片，视频等均是来源于用户投稿和互联网及文摘转载整编而成，不代表本站观点，不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益，请在线联系站长,一经查实,本站将立刻删除。本文来自网络,若有侵权，请联系删除，如若转载，请注明出处：https://haidsoft.com/113631.html

Python库：gdal