地理空间数据格式转化库之GDAL

随笔 2017-02-12

  GDAL是一个用于光栅和矢量地理空间数据格式转换的开源库。如果要读取空间数据这个库是必不可少的,它可以操作各种栅格地理数据格式,包括读取、写入、转换、处理各种栅格数据。这个库还同时包括了操作矢量数据的另一个有名的库OGR,这样这个库就同时具备了操作栅格矢量数据的能力。

  GDAL支持跨平台和多语言,C++版的编译配置网上有很多资料,但也很容易遇到错误,因此推荐使用Python版。使用pip install gdal安装即可,如果用Anconda科学计算包也可以使用conda install gdal安装。

  1.起步:读取并输出图像信息

 from osgeo import gdal

file_path = 'image.tif' //change to your own path
dataset = gdal.Open(file_path)
print(type(dataset))

metadata = dataset.GetMetadata()
print(metadata)

projection = dataset.GetProjection()
print(projection)

  2.基础:几何图像操作

from osgeo import ogr  
multiline = ogr.Geometry(ogr.wkbMultiLineString)

line1 = ogr.Geometry(ogr.wkbLineString)
line1.AddPoint(1214242.4174581182, 617041.9717021306)
line1.AddPoint(1234593.142744733, 629529.9167643716)
multiline.AddGeometry(line1)

line1 = ogr.Geometry(ogr.wkbLineString)
line1.AddPoint(1184641.3624957693, 626754.8178616514)
line1.AddPoint(1219792.6152635587, 606866.6090588232)
multiline.AddGeometry(line1)

print multiline.ExportToWkt()

  简单的图形只需要构建并添加点即可。

  3.矢量图层操作

  主要用到OGR模块,对常用的.shp,.osm格式文件进行操作,即地理数据常用的格式。 同时支持从PostgreSQL和网络获取数据。

import os
from osgeo import ogr

daShapefile = r"C:\Temp\Voting_Centers_and_Ballot_Sites.shp"

driver = ogr.GetDriverByName('ESRI Shapefile')
#需要用到驱动器
dataSource = driver.Open(daShapefile, 0) # 0 表示只读. 1 表示可写.

# 检测文件是否存在
if dataSource is None:
    print 'Could not open %s' % (daShapefile)
else:
    print 'Opened %s' % (daShapefile)
    layer = dataSource.GetLayer()
    featureCount = layer.GetFeatureCount()
    print "Number of features in %s: %d" % (os.path.basename(daShapefile),featureCount)

  4.光栅图层操作

  支持将OGR数据转化为光栅数据

   from osgeo import gdal, ogr

# 定义新栅格的像素和NoData值
pixel_size = 25
NoData_value = -9999

# 输入文件
vector_fn = 'test.shp'

# 新的tif文件被创建
raster_fn = 'test.tif'

# 打开数据源
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 创建目标数据源
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

#栅格化
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])

  以上只是GDAL的一小部分功能,它对于地理数据的强大操作能力还有待进一步探索。

附:Python版的文档


本文由 Tony 创作,采用 知识共享署名 3.0,可自由转载、引用,但需署名作者且注明文章出处。

赏个馒头吧