GDAL工具集介绍

GDAL是一个库,但同时它的安装程序里又带有很多实用的工具程序(这些程序在Debian中被称为gdal-bin)。 在此把这些工具列一下,因为在后面的章节中会用到这些工具。 但是因为有些概念用到了后面的内容,所以在此不进行详细的说明。

GDAL附带工具

GDAL创建了下面的工具程序

  • gdal2tiles.py
  • gdal2xyz.py
  • gdaladdo - 为文件添加影像金字塔。
  • gdalbuildvrt
  • gdalchksum.py
  • gdal-config - 获取建立使用gdal的软件需要的参数。
  • gdal_contour - 给DEM创建等高线。
  • gdalenhance - 进行图像增强。
  • gdal_grid
  • gdalident.py
  • gdalimport.py
  • gdalinfo - 统计文件的信息。
  • gdalmanage
  • gdal_merge.py - 把一系列的图像创建一个快速的地图集。
  • gdal_polygonize.py
  • gdal_proximity.py
  • gdal_rasterize - 把一个矢量图层栅格化。
  • gdal_retile.py
  • gdal_sieve.py
  • gdaltindex - 为MapServer建立一个栅格索引。
  • gdaltransform
  • gdal_translate - 拷贝一个可以控制输出格式的栅格文件。
  • gdalwarp - 把一个图像转换到一个新的坐标系统。
  • nearblack - 转换临近的黑/白边框为精准数值。
  • pct2rgb.py - 把一个8位调色板图像转化成24位位图。
  • rgb2pct.py - 把一个24位的图像转化成8位调色板图像。

总的命令行参数

所有的GDAL命令行工具程序都支持下面的“总的”命令行参数。 这些参数列出的是GDAL工具的信息,而并不是某个具体GDAL程序的信息。

  • --version 登记GDAL版本并退出。
  • --formats 列出所有GDAL支持的栅格格式(只读和读写的)并退出。 ro是只读驱动,rw是读写驱动(如支持CreateCopy()); rw+是读写和更新驱动(如支持Create())。
  • --format format 列出单个格式驱动的细节信息。 格式名需要是在--formats 中列出的格式名。比如GTiff。
  • --optfile file 读取指定名称的文件并把其中的内容当成参数传入命令行列表。 行首以#开头的行将被忽略。
  • --config key value 设置配置把指定键设置为某个值,而不必把他们设置为环境变量。 一些命令参数键是GDAL_CACHEMAX(用于缓存的内存有多少M) 以及GDAL_DATA(gdal的数据路径)。单一的驱动会被其他配置参数影响。
  • --debug value 控制调试信息的打印。ON值表示允许调试信息输出, OFF值表示不要输出调试信息。
  • --help-general 给出一个简短的普通的GDAL命令行参数用法信息然后退出。
$ gdalinfo --version
GDAL 1.6.3, released 2009/11/19
$ gdal_translate --version
GDAL 1.6.3, released 2009/11/19
$ gdal_translate --formats
Supported Formats:
  GRASS (ro): GRASS Database Rasters (5.7+)
  VRT (rw+): Virtual Raster
  GTiff (rw+): GeoTIFF
......

创建新的文件

存取一个已存在的文件来读取是一件很容易的事情, 只需要在命令行中指定文件或者数据集的名字。 但是,创建一个文件是一件非常复杂的事情。 你可能需要指定创建格式,各种创建参数,以及指定一个坐标系统。 在不同的GDAL工具中有许多参数都是差不多的,列举一下:

  • -of format 选择要创建新的文件的格式。 这个格式被指定为类似GTiff(GeoTIFF格式)或者HFA(ERDAS格式),格式的定义与前面是一致的。 所有的格式列表可以用--formats 参数列出来。 只有格式列表rw可以被写入。默认是创建GeoTIFF文件。
  • -co NAME=VALUE 许多格式会有一个或者更多的创建参数来控制文件创建的细节。 比如GeoTIFF可以用创建参数控制压缩, 或者控制是否用分片还是分带来进行存储。 可以使用的创建参数根据格式驱动不同而不同, 而一些简单的格式根本就没有创建参数。
  • -a_srs SRS 有几个工具(gdal_translategdalwarp)可以在命令行中 通过类似-a_srs(分配输出SRS)、 -s_srs(源SRS)、-t_srs(目标SRS)来指定坐标系统。 这些工具允许以一系列格式定义坐标系统(SRS就是空间参考系统), 关于坐标系统的介绍,见第[chp:proj]章。

gdalinfo

列出栅格数据集的信息

gdalinfo程序输出gdal支持的栅格格式的一系列信息:

  • 用于存取的文件格式驱动
  • 栅格大小(行列数)
  • 文件的坐标系统(OGC WKT形式)
  • 图像关联到地理的转换参数(当前不包含旋转系数)
  • 地理上的边界坐标,如果可能的话还有基于经纬度的完整的地理转换参数(如果是GCPs就没有)
  • 地面控制点(GCPs)
  • 广义的(包括子栅格的原数据)文件原数据
  • 波段数据类型
  • 波段光度解析
  • 波段瓦片大小
  • 波段描述
  • 波段最大最小值(已经经过计算的)
  • 波段无意义值
  • 波段可获得的略缩图分辨率
  • 波段单位类型(如:波段的高程是米制还是英制)
  • 波段的假颜色列表

例子:

$ gdalinfo ~/openev/utm.tif

这个输出就不在此列出了。

gdal_translate

gdal_translate 工具可以用来在不同格式间转换栅格数据。 并且在处理过程中做一些诸如子栅格设置、重采样和像元值集体变化等工作。

例子:

$ gdal_translate -of GTiff -co "TILED=YES" utm.tif utm_tiled.tif

gdalenhance

进行图像增强。增强类型通过使用 -equalize,或其他参数来实现。 -s_nodata参数设置图像中认为是空值的数值。

$ gdalenhance -s_nodata 0 -equalize

gdaladdo

gdaladdo用于建立或者重建图像金字塔。

对于没有图像金字塔的文件,使用ArcCatalog或ArcMap进行浏览时, 会要求生成金字塔文件。然后在同一目录下生成文件名相同的rrd文件。

gdaladdo工具可以为大多数支持的格式建立或者重建金字塔。 可以使用下面的重采样算法来进行缩小重采样操作。

例如选择一个缩放水平如2表示略缩图缩放程度是源图像每个维上分辨率的1/2的缩略图。 若文件在所选缩放水平上已经存在略缩图,那么那个缩放程度将被重新计算并覆盖写入。

一些格式根本不支持金字塔。许多格式在文件以外以扩展名.ovr存储金字塔,TIFF就是如此。 GeoTIFF格式直接把金字塔存储到原有的文件中。

在TIFF中创建金字塔可以通过用COMPRESS_OVERVIEW配置参数进行压缩。 所有GeoTIFF支持的压缩方法,可以在这里获得(如: --config COMPRESS_OVERVIEW DEFLATE)。

大多数驱动也支持一个备用的略缩图格式(使用的是Erdas图像格式)。 引发这个备用格式使用 USE_RRD=YES 来设置参数。 这样做会把GDAL程序创建的金字塔放到一个辅助的.aux文件中使得可以该金字塔可以直接在Erdas中使用或者也可以在ArcGIS中使用。

例子:

在所提供的TIFF文件内部创建金字塔

$ gdaladdo -r average abc.tif 2 4 8 16

从一个ERDAS.IMG文件中创建一个外部的压缩的金字塔文件。

$ gdaladdo --config COMPRESS\_OVERVIEW DEFLATE erdas.img 2 4 8 16

为给定JPEG文件创建一个Erdas Imagine 格式金字塔

$ gdaladdo --config USE\_RRD YES airphoto.jpg 3 9 27 81

gdal_merge

[style=MyBash]

$ gdal_merge -o d:merge150421.tif -ps 30 30 -n 0 k50012.tif k50024.tif k50023.tif

要注意这个图像的叠放顺序为“后来者居上”。

gdalwarp

gdalwarp工具是一个图像镶嵌、重投影、和绑定的工具,此工具默认输出为BIP存储格式的数据。 程序可以重投影到任何支持的投影,而且如果图像(“raw” with)控制信息也可以把GCPs和图像存储在一起。

例子:

一个带有用经纬度标记的边界控制点的8位GeoTiff格式的Spot影像可以通过下面的命令绑定到一个UTM投影上。

$ gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif

一个带有用经纬度标记的边界控制点的两个波段的HDF格式的ASTER影像可以通过下面的命令绑定到UTM投影上。

$ gdalwarp HDF4_SDS:ASTER_L1B:"pg-PR1B0000-2002031402_100_001":2 pg-PR1B0000-2002031402_100_001_2.tif

下面假定为中国常用高斯——克吕格投影,投影文件定义为PROJ.4格式。以51带为例:

+proj=tmerc
+lat_0=0
+lon_0=123
+k=1
+x_0=21500000
+y_0=0
+ellps=krass

其他的50带,52带,53带请自行定义。将其投影转换为Albers投影:

+proj=aea
+lat_1=25
+lat_2=47
+lat_0=0
+lon_0=105
+x_0=0
+y_0=0
+ellps=krass

这两个文件分别所设为g51.pj4与al.pj4。那么转换命令为:

$ gdalwarp -s_srs g51.pj4 -t_srs al.pj4 l51036p_a.tif wl50036p_a.tif

栅格数据的投影

首先要知道输入投影和输出投影的WKT(Well Known Text),可以通过 GetProjection() 读到, 可以用 SpatialReference 对象创建。

用下面的语句新建栅格数据集并重新投影,投影结果输出到新数据集之中:

gdal.CreateAndReprojectImage( <source_dataset>, <output_filename>, src_wkt=<source_wkt>,       dst_wkt=<output_wkt>, dst_driver=<Driver>, eResampleAlg = <GDALResampleAlg>)

对比两种方法的执行效率:逐个像素处理,使用内建函数处理, 内建的函数效率最高,比逐个像素循环要快很多。

GDAL和 Pillow 的互操作

前面提到过,GDAL和PIL很相像。它们处理和操作的对象都是栅格图像。 但它们又不一样。 GDAL主要重点放在地理或遥感数据的读写和数据建模以及地理定位和转换, 但是PIL的重点是放在图像本身处理上的。

至于在底层数据处理上,两者都可以用 numpy 转化的二进制作为数据处理。 所以,理论上是可以互相共享和交换数据的。实际上也确实可以。

首先,我要说明的是GDAL的核心在波段(band), 一切操作的基础和核心都在波段。 波段可以单独拿出来操作,至于波段在数据集中的顺序无关紧要。 因为遥感图像大多比RGB图像的波段要多,而每个波段单独都是一个完整的整体, 每个波段单独拿出来都是一个数据集。而 Pillow 的核心在数据集(DataSet)这里的概念是对应GDAL中的数据集的概念。 当然,在 Pillow 本身中没有这种说法,也就是不把波段单独操作, 操作大部分需要RGB一体化地进行。

两部分的操作的主要衔接部分就是创建、读取与写入。 读取数据后怎么处理是两个库各自的事情。 所以这里主要内容就是介绍两个库各自的创建,读取和写入的操作,以及两个库的过渡。

使用GDAL读取数据

比较两个库的读取,GDAL读取一个图像中的数据

In [15]:
from osgeo import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
data_arr = dataset.ReadAsArray(30,70,5,5)
type(data_arr)
Out[15]:
numpy.ndarray
In [16]:
data_arr
Out[16]:
array([[[147, 141, 151, 146, 145],
        [148, 149, 151, 143, 139],
        [163, 164, 162, 152, 149],
        [167, 169, 164, 160, 159],
        [168, 172, 162, 162, 164]],

       [[  7,   4,  17,  12,  11],
        [  7,  10,  14,   6,   2],
        [ 10,  11,  11,   3,   0],
        [  8,  10,   8,   4,   4],
        [ 12,  16,   6,   6,   9]],

       [[ 18,  12,  24,  19,  18],
        [ 16,  17,  21,  13,   9],
        [ 15,  16,  16,   7,   6],
        [ 13,  14,  11,   8,  10],
        [ 16,  20,  10,  10,  15]]], dtype=uint8)
In [17]:
data_bin = dataset.ReadRaster(30,70,5,5)
data_bin
Out[17]:
b'\x93\x8d\x97\x92\x91\x94\x95\x97\x8f\x8b\xa3\xa4\xa2\x98\x95\xa7\xa9\xa4\xa0\x9f\xa8\xac\xa2\xa2\xa4\x07\x04\x11\x0c\x0b\x07\n\x0e\x06\x02\n\x0b\x0b\x03\x00\x08\n\x08\x04\x04\x0c\x10\x06\x06\t\x12\x0c\x18\x13\x12\x10\x11\x15\r\t\x0f\x10\x10\x07\x06\r\x0e\x0b\x08\n\x10\x14\n\n\x0f'

打开了数据集,就有两种方法来获取数据。

虽然读出的一个是二进制,一个是数组,Numeric数组用tostirng转换出来的二进制和用ReadAsArray读出的相同。

In [18]:
data_arr.tostring()
Out[18]:
b'\x93\x8d\x97\x92\x91\x94\x95\x97\x8f\x8b\xa3\xa4\xa2\x98\x95\xa7\xa9\xa4\xa0\x9f\xa8\xac\xa2\xa2\xa4\x07\x04\x11\x0c\x0b\x07\n\x0e\x06\x02\n\x0b\x0b\x03\x00\x08\n\x08\x04\x04\x0c\x10\x06\x06\t\x12\x0c\x18\x13\x12\x10\x11\x15\r\t\x0f\x10\x10\x07\x06\r\x0e\x0b\x08\n\x10\x14\n\n\x0f'

从波段中获取数据和从数据集中获取数据的方法十分相似。

使用Pillow读取数据

注意,使用Pillow读取时,要注意其类型。

from PIL import Image im = Image.open(’K52E015007.tif’) region = im.crop((30,70,35,75)) region.tostring()

In [19]:
from PIL import Image
im = Image.open("/gdata/geotiff_file.tif")
In [20]:
region = im.crop((30,70,35,75))
region.tobytes()
Out[20]:
b'\x93\x07\x12\x8d\x04\x0c\x97\x11\x18\x92\x0c\x13\x91\x0b\x12\x94\x07\x10\x95\n\x11\x97\x0e\x15\x8f\x06\r\x8b\x02\t\xa3\n\x0f\xa4\x0b\x10\xa2\x0b\x10\x98\x03\x07\x95\x00\x06\xa7\x08\r\xa9\n\x0e\xa4\x08\x0b\xa0\x04\x08\x9f\x04\n\xa8\x0c\x10\xac\x10\x14\xa2\x06\n\xa2\x06\n\xa4\t\x0f'

im可以类比成gdal的dataset,im也可以从DataSet中提取某个范围的数据。

可以看出,虽然读取的都是同样位置的数据,但是输出的结果不一样。

Pillow与GDAL读取数据的转换

这里注意,GDAL与Pillow的空间模型并不一致。 在Pillow的下,截取区域矩形的定义和GDAL不同,GDAL是顶点X、顶点Y、宽、高; Pillow是顶点X、顶点Y、终点X,终点Y。 这就是GDAL和Pillow的区别。转换一下:

import numpy as np data = dataset.ReadAsArray(30,70,5,5) datas = [i for i in data] from numpy import reshape datas = [reshape(i,(-1,1)) for i in data] datas = np.concatenate(datas,1) datas.tostring()

In [21]:
import numpy as np
data = dataset.ReadAsArray(30,70,5,5)
datas = [i for i in data]
from numpy import reshape
datas = [reshape(i,(-1,1)) for i in data]
datas = np.concatenate(datas,1)
datas.tostring()
Out[21]:
b'\x93\x07\x12\x8d\x04\x0c\x97\x11\x18\x92\x0c\x13\x91\x0b\x12\x94\x07\x10\x95\n\x11\x97\x0e\x15\x8f\x06\r\x8b\x02\t\xa3\n\x0f\xa4\x0b\x10\xa2\x0b\x10\x98\x03\x07\x95\x00\x06\xa7\x08\r\xa9\n\x0e\xa4\x08\x0b\xa0\x04\x08\x9f\x04\n\xa8\x0c\x10\xac\x10\x14\xa2\x06\n\xa2\x06\n\xa4\t\x0f'

可以看到现在结果一致了。

这里就表现了两个库的设计概念模型的不同。 GDAL把图像看成是由不同传感器获取的不同频率的电磁波构成的影像文件,读取的数据是默认的以band组织的; Pillow则把图像看成是由单个像素构成的,每个像素是记录的由RGB三色构成的像素颜色的数据。

从波段来看

如果是单个波段,就不存在RGB存储的问题了。使用下面的方式打开,可以看出读取数据时,两个库读取的结果是一样的。

r,g,b = region.split() r.tostring() band = dataset.GetRasterBand(1) band.ReadRaster(30,70,5,5)

In [22]:
r,g,b = region.split()
r.tobytes()  
Out[22]:
b'\x93\x8d\x97\x92\x91\x94\x95\x97\x8f\x8b\xa3\xa4\xa2\x98\x95\xa7\xa9\xa4\xa0\x9f\xa8\xac\xa2\xa2\xa4'
In [23]:
band = dataset.GetRasterBand(1)
band.ReadRaster(30,70,5,5)
Out[23]:
b'\x93\x8d\x97\x92\x91\x94\x95\x97\x8f\x8b\xa3\xa4\xa2\x98\x95\xa7\xa9\xa4\xa0\x9f\xa8\xac\xa2\xa2\xa4'

从写数据来看

再比较两个库的写入,写入数据gdal用的是WriteRaster。同样DataSet一个,Band一个。

help(ds.WriteRaster)
help(band.WriteRaster)
help(band.WriteArray)

这个就是把一个二进制字符串放到buf_string位置, 那个xoff,yoff,xsize,ysize就是要把数据贴到整张图的什么位置, buf_xsize,buf_ysize是那个二进制字符串表示的区域大小 (建议和xsize,ysize一致,不然gdal会自动缩放, 而且自动缩放的方法是我们不能控制的,只会用最临近法)。 buf_type说明的是输入的二进制字符串的数据类型, 如果和原底图的数据类型不一样,会自动扩充和截断数据长度, band_list就是要输入的波段的列表。

PIL中对数据的写入用的是paste,

help(im.paste)

好了。既然上面验证了两个库中读取的数据二进制一样,那么就可以互相交换数据了。 可以把gdal的数据读出,贴到pil打开的数据中, 也可以把PIL的数据读出,贴到gdal打开的数据中。

在PIL中还有一个很实用的代码--fromstring

help(Image.fromstring)

Image本身有一个,im(Image打开的对象)也有一个。 创建一个空数据,然后用个fromstring, 然后保存后就生成一张新图! 保存一下刚才那个10*10的图,看看是什么……

newim = Image.new("RGB",(10,10))
newim.fromstring(datas.tostring())
newim.save("f:/png/new.png","png")

另一种:

Image.fromstring("RGB",(10,10),datas.tostring()).save("f:/png/new1.png","png")

只要一行!