gdal和PIL的互操作
前面提到过，gdal和PIL很相像。它们处理和操作的对象都是栅格图像。但它们又不一样。gdal主要重点放在地理或遥感数据的读写和数据建模以及地理定位和转换，但是PIL的重点是放在图像本身处理上的。
至于在底层数据处理上，两者都可以用Numeric转化的二进制作为数据处理。所以，理论上是可以互相共享和交换数据的。实际上也确实可以。
首先，我要说明的是gdal的核心在波段（band），一切操作的基础和核心都在波段。波段可以单独拿出来操作，至于波段在数据集中的顺序无关紧要。因为遥感图像大多比RGB图像的波段要多，而每个波段单独都是一个完整的整体，每个波段单独拿出来就是一个数据集。而PIL的核心在数据集（DataSet）这里数据集的概念是对应gdal中的数据集的概念，可能在PIL本身中没有这种说法。也就是不把波段单独操作。操作大部分需要RGB一体化地进行，虽然可以单波段操作，但是我想一般人很少这样干！
两部分的操作的主要衔接部分就是创建和读取，以及写入。至于进去了里面怎么变化就是两个库各自的事情了。所以这篇文章的主要内容就是介绍两个库各自的创建，读取和写入的操作，以及两个库的过渡。
比较两个库的读取，gdal读取一个图像中的数据
打开了数据集，就有两种方法来获取数据。一种是通过数据集
1 >>> help(dataset.ReadRaster)
2 Help on method ReadRaster in module gdal:
3 ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t
4 ype=None, band_list=None) method of gdal.Dataset instance
5 >>> help(dataset.ReadAsArray)
6 Help on method ReadAsArray in module gdal:
7 ReadAsArray(self, xoff=0, yoff=0, xsize=None, ysize=None) method of gdal.Dataset
8 instance
9 >>>
其中ReadRaster读出的是二进制（在python是str类型），ReadAsArray读出的是Numeric数组。
虽然读出的一个是二进制，一个是数组，Numeric数组用tostirng转换出来的二进制和用ReadAsArray读出的相同。
>>> ds.ReadAsArray(130,90,10,10).tostring()
'\xd4\xd8\xd4\xcf\xcc\xce\xce\xc3\x7f\xbf\xd1\xd5\xca\xd4\xd3\xd3\xd1\xc0|\xc5\x
cd\xcd\xd1\xd1\xca\xd0\xce\xbf\x88\xc3\xd0\xd0\xce\xcb\xca\xca\xcc\xb8w\xb9\xcc\
xcf\xca\xca\xd4\xd4\xcb\xb5\x87\xbe\xd3\xd4\xd0\xcd\xce\xbf\xcf\xc6\x80\xc1\xd3\
xd3\xd0\xcb\xcc\xcc\xd0\xc0\x86\xbe\xd0\xce\xd1\xce\xd5\xce\xd0\xc4\x88\xc5\xd0\
xce\xc9\xc9\xd4\xd1\xd0\xc3\x85\xd0\xce\xd2\xc9\xcc\xcc\xd2\xd0\xcd\xa8\xcb\xce\
xd2\xcf\xca\xc7\xcb\xce\xc5\x83\xc1\xcb\xcf\xc7\xd1\xd0\xd1\xd1\xc2\x80\xc7\xc7\
xc7\xcf\xcf\xc8\xce\xce\xc1\x8c\xc5\xca\xca\xcc\xc8\xc8\xc8\xcc\xba{\xbb\xc8\xc9
\xc7\xc2\xd0\xd2\xc9\xb5\x88\xc0\xcf\xce\xcd\xc5\xca\xbd\xcc\xc6\x82\xc3\xd0\xcd
\xcd\xc3\xc8\xc9\xcd\xc0\x88\xc0\xca\xc8\xcb\xc8\xcf\xca\xcc\xc0\x84\xc1\xca\xc8
\xc5\xc5\xd0\xcd\xcc\xbf\x81\xcc\xc8\xcc\xc7\xca\xca\xcf\xcd\xc9\xa4\xc8\xd3\xd7
\xcd\xc8\xc5\xcd\xd3\xcd\x8d\xc9\xd0\xd4\xc7\xd1\xd0\xd5\xd6\xca\x8a\xcf\xcc\xcc
\xd1\xd1\xca\xd2\xd3\xc9\x96\xcd\xcf\xcf\xce\xca\xca\xcc\xd1\xc2\x85\xc3\xcd\xce
\xc9\xc8\xd6\xd6\xce\xbe\x93\xc6\xd4\xd3\xcf\xcb\xd0\xc2\xd4\xcd\x8a\xc9\xd2\xd0
\xcf\xc9\xce\xd3\xd5\xc7\x8c\xc6\xcb\xc9\xd0\xcd\xd4\xcf\xd2\xc8\x8c\xc7\xcb\xc9
\xca\xca\xd5\xd2\xd1\xc5\x89\xd1\xc9\xcd\xcb\xce\xce\xd1\xcf\xce\xaa\xca'
>>> ds.ReadRaster(130,90,10,10)
'\xd4\xd8\xd4\xcf\xcc\xce\xce\xc3\x7f\xbf\xd1\xd5\xca\xd4\xd3\xd3\xd1\xc0|\xc5\x
cd\xcd\xd1\xd1\xca\xd0\xce\xbf\x88\xc3\xd0\xd0\xce\xcb\xca\xca\xcc\xb8w\xb9\xcc\
xcf\xca\xca\xd4\xd4\xcb\xb5\x87\xbe\xd3\xd4\xd0\xcd\xce\xbf\xcf\xc6\x80\xc1\xd3\
xd3\xd0\xcb\xcc\xcc\xd0\xc0\x86\xbe\xd0\xce\xd1\xce\xd5\xce\xd0\xc4\x88\xc5\xd0\
xce\xc9\xc9\xd4\xd1\xd0\xc3\x85\xd0\xce\xd2\xc9\xcc\xcc\xd2\xd0\xcd\xa8\xcb\xce\
xd2\xcf\xca\xc7\xcb\xce\xc5\x83\xc1\xcb\xcf\xc7\xd1\xd0\xd1\xd1\xc2\x80\xc7\xc7\
xc7\xcf\xcf\xc8\xce\xce\xc1\x8c\xc5\xca\xca\xcc\xc8\xc8\xc8\xcc\xba{\xbb\xc8\xc9
\xc7\xc2\xd0\xd2\xc9\xb5\x88\xc0\xcf\xce\xcd\xc5\xca\xbd\xcc\xc6\x82\xc3\xd0\xcd
\xcd\xc3\xc8\xc9\xcd\xc0\x88\xc0\xca\xc8\xcb\xc8\xcf\xca\xcc\xc0\x84\xc1\xca\xc8
\xc5\xc5\xd0\xcd\xcc\xbf\x81\xcc\xc8\xcc\xc7\xca\xca\xcf\xcd\xc9\xa4\xc8\xd3\xd7
\xcd\xc8\xc5\xcd\xd3\xcd\x8d\xc9\xd0\xd4\xc7\xd1\xd0\xd5\xd6\xca\x8a\xcf\xcc\xcc
\xd1\xd1\xca\xd2\xd3\xc9\x96\xcd\xcf\xcf\xce\xca\xca\xcc\xd1\xc2\x85\xc3\xcd\xce
\xc9\xc8\xd6\xd6\xce\xbe\x93\xc6\xd4\xd3\xcf\xcb\xd0\xc2\xd4\xcd\x8a\xc9\xd2\xd0
\xcf\xc9\xce\xd3\xd5\xc7\x8c\xc6\xcb\xc9\xd0\xcd\xd4\xcf\xd2\xc8\x8c\xc7\xcb\xc9
\xca\xca\xd5\xd2\xd1\xc5\x89\xd1\xc9\xcd\xcb\xce\xce\xd1\xcf\xce\xaa\xca'
>>>从波段中获取数据和从数据集中获取数据有十分相似的方法。
1 >>> help(band.ReadAsArray)
2 Help on method ReadAsArray in module gdal:
3 ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None
4 , buf_ysize=None, buf_obj=None) method of gdal.Band instance
5 >>> help(band.ReadRaster)
6 Help on method ReadRaster in module gdal:
7 ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t
8 ype=None) method of gdal.Band instance
9 >>>
而PIL打开获取数据的过程如下：
1 >>> import Image
2 >>> im = Image.open("E:/gisdata/gtif/sd.tif")
3 >>> dir(im)
4 ['_Image__transformer', '_TiffImageFile__first', '_TiffImageFile__fp', '_TiffIma
5 geFile__frame', '_TiffImageFile__next', '__doc__', '__init__', '__module__', '_c
6 ompression', '_copy', '_decoder', '_dump', '_expand', '_makeself', '_new', '_ope
7 n', '_planar_configuration', '_seek', '_setup', '_tell', 'category', 'convert',
8 'copy', 'crop', 'decoderconfig', 'decodermaxblock', 'draft', 'filename', 'filter
9 ', 'format', 'format_description', 'fp', 'fromstring', 'getbands', 'getbbox', 'g
10 etcolors', 'getdata', 'getextrema', 'getim', 'getpalette', 'getpixel', 'getproje
11 ction', 'histogram', 'ifd', 'im', 'info', 'load', 'load_end', 'load_prepare', 'm
12 ode', 'offset', 'palette', 'paste', 'point', 'putalpha', 'putdata', 'putpalette'
13 , 'putpixel', 'quantize', 'readonly', 'resize', 'rotate', 'save', 'seek', 'show'
14 , 'size', 'split', 'tag', 'tell', 'thumbnail', 'tile', 'tobitmap', 'tostring', '
15 transform', 'transpose', 'verify']
im可以类比成gdal的dataset，im也可以从DataSet中提取某个范围的数据。
>>> region = im.crop((130,90,140,100))
>>> region.tostring()
'\xd4\xce\xd3\xd8\xd2\xd7\xd4\xcf\xcd\xcf\xca\xc8\xcc\xc7\xc5\xce\xcb\xcd\xce\xc
e\xd3\xc3\xc5\xcd\x7f\x83\x8d\xbf\xc1\xc9\xd1\xcb\xd0\xd5\xcf\xd4\xca\xc7\xc7\xd
4\xd1\xd1\xd3\xd0\xd0\xd3\xd1\xd5\xd1\xd1\xd6\xc0\xc2\xca|\x80\x8a\xc5\xc7\xcf\x
cd\xc7\xcc\xcd\xc7\xcc\xd1\xcf\xd1\xd1\xcf\xd1\xca\xc8\xca\xd0\xce\xd2\xce\xce\x
d3\xbf\xc1\xc9\x88\x8c\x96\xc3\xc5\xcd\xd0\xca\xcf\xd0\xca\xcf\xce\xcc\xce\xcb\x
c8\xca\xca\xc8\xca\xca\xc8\xcc\xcc\xcc\xd1\xb8\xba\xc2w{\x85\xb9\xbb\xc3\xcc\xc8
\xcd\xcf\xc9\xce\xca\xc7\xc9\xca\xc2\xc8\xd4\xd0\xd6\xd4\xd2\xd6\xcb\xc9\xce\xb5
\xb5\xbe\x87\x88\x93\xbe\xc0\xc6\xd3\xcf\xd4\xd4\xce\xd3\xd0\xcd\xcf\xcd\xc5\xcb
\xce\xca\xd0\xbf\xbd\xc2\xcf\xcc\xd4\xc6\xc6\xcd\x80\x82\x8a\xc1\xc3\xc9\xd3\xd0
\xd2\xd3\xcd\xd0\xd0\xcd\xcf\xcb\xc3\xc9\xcc\xc8\xce\xcc\xc9\xd3\xd0\xcd\xd5\xc0
\xc0\xc7\x86\x88\x8c\xbe\xc0\xc6\xd0\xca\xcb\xce\xc8\xc9\xd1\xcb\xd0\xce\xc8\xcd
\xd5\xcf\xd4\xce\xca\xcf\xd0\xcc\xd2\xc4\xc0\xc8\x88\x84\x8c\xc5\xc1\xc7\xd0\xca
\xcb\xce\xc8\xc9\xc9\xc5\xca\xc9\xc5\xca\xd4\xd0\xd5\xd1\xcd\xd2\xd0\xcc\xd1\xc3
\xbf\xc5\x85\x81\x89\xd0\xcc\xd1\xce\xc8\xc9\xd2\xcc\xcd\xc9\xc7\xcb\xcc\xca\xce
\xcc\xca\xce\xd2\xcf\xd1\xd0\xcd\xcf\xcd\xc9\xce\xa8\xa4\xaa\xcb\xc8\xca'
>>>这里注意，在pil的下，截取区域矩形的定义和gdal不同，gdal是顶点X、顶点Y、宽、高；pil是顶点X、顶点Y、终点X，终点Y
诶，似乎不一样啊！
确实不一样，不过这就是gdal和pil的区别了。转换一下吧。
激动人心的时刻到了！当当当当～～～～～
>>> datas.tostring()
'\xd4\xce\xd3\xd8\xd2\xd7\xd4\xcf\xcd\xcf\xca\xc8\xcc\xc7\xc5\xce\xcb\xcd\xce\xc
e\xd3\xc3\xc5\xcd\x7f\x83\x8d\xbf\xc1\xc9\xd1\xcb\xd0\xd5\xcf\xd4\xca\xc7\xc7\xd
4\xd1\xd1\xd3\xd0\xd0\xd3\xd1\xd5\xd1\xd1\xd6\xc0\xc2\xca|\x80\x8a\xc5\xc7\xcf\x
cd\xc7\xcc\xcd\xc7\xcc\xd1\xcf\xd1\xd1\xcf\xd1\xca\xc8\xca\xd0\xce\xd2\xce\xce\x
d3\xbf\xc1\xc9\x88\x8c\x96\xc3\xc5\xcd\xd0\xca\xcf\xd0\xca\xcf\xce\xcc\xce\xcb\x
c8\xca\xca\xc8\xca\xca\xc8\xcc\xcc\xcc\xd1\xb8\xba\xc2w{\x85\xb9\xbb\xc3\xcc\xc8
\xcd\xcf\xc9\xce\xca\xc7\xc9\xca\xc2\xc8\xd4\xd0\xd6\xd4\xd2\xd6\xcb\xc9\xce\xb5
\xb5\xbe\x87\x88\x93\xbe\xc0\xc6\xd3\xcf\xd4\xd4\xce\xd3\xd0\xcd\xcf\xcd\xc5\xcb
\xce\xca\xd0\xbf\xbd\xc2\xcf\xcc\xd4\xc6\xc6\xcd\x80\x82\x8a\xc1\xc3\xc9\xd3\xd0
\xd2\xd3\xcd\xd0\xd0\xcd\xcf\xcb\xc3\xc9\xcc\xc8\xce\xcc\xc9\xd3\xd0\xcd\xd5\xc0
\xc0\xc7\x86\x88\x8c\xbe\xc0\xc6\xd0\xca\xcb\xce\xc8\xc9\xd1\xcb\xd0\xce\xc8\xcd
\xd5\xcf\xd4\xce\xca\xcf\xd0\xcc\xd2\xc4\xc0\xc8\x88\x84\x8c\xc5\xc1\xc7\xd0\xca
\xcb\xce\xc8\xc9\xc9\xc5\xca\xc9\xc5\xca\xd4\xd0\xd5\xd1\xcd\xd2\xd0\xcc\xd1\xc3
\xbf\xc5\x85\x81\x89\xd0\xcc\xd1\xce\xc8\xc9\xd2\xcc\xcd\xc9\xc7\xcb\xcc\xca\xce
\xcc\xca\xce\xd2\xcf\xd1\xd0\xcd\xcf\xcd\xc9\xce\xa8\xa4\xaa\xcb\xc8\xca'
>>>相同了吧！
这里就表现了两个库的设计哲学不同了。gdal的核心是band，读取的数据是默认以band组织的，pil的核心是dataset，读取的数据默认是统一组织的。换句话说，gdal看比较顺眼，pil看{RGB,RGB,RGB……}比较顺眼。
同时把数据降低到Band水平（其实不能说降低，Band的水平不比DataSet低），就可以看出来在Band水平，两个库是一样一样一样的啊！
>>> r,g,b = region.split() >>> r.tostring() '\xd4\xd8\xd4\xcf\xcc\xce\xce\xc3\x7f\xbf\xd1\xd5\xca\xd4\xd3\xd3\xd1\xc0|\xc5\x cd\xcd\xd1\xd1\xca\xd0\xce\xbf\x88\xc3\xd0\xd0\xce\xcb\xca\xca\xcc\xb8w\xb9\xcc\ xcf\xca\xca\xd4\xd4\xcb\xb5\x87\xbe\xd3\xd4\xd0\xcd\xce\xbf\xcf\xc6\x80\xc1\xd3\ xd3\xd0\xcb\xcc\xcc\xd0\xc0\x86\xbe\xd0\xce\xd1\xce\xd5\xce\xd0\xc4\x88\xc5\xd0\ xce\xc9\xc9\xd4\xd1\xd0\xc3\x85\xd0\xce\xd2\xc9\xcc\xcc\xd2\xd0\xcd\xa8\xcb' >>> band = ds.GetRasterBand(1) >>> band.ReadRaster(130,90,10,10) '\xd4\xd8\xd4\xcf\xcc\xce\xce\xc3\x7f\xbf\xd1\xd5\xca\xd4\xd3\xd3\xd1\xc0|\xc5\x cd\xcd\xd1\xd1\xca\xd0\xce\xbf\x88\xc3\xd0\xd0\xce\xcb\xca\xca\xcc\xb8w\xb9\xcc\ xcf\xca\xca\xd4\xd4\xcb\xb5\x87\xbe\xd3\xd4\xd0\xcd\xce\xbf\xcf\xc6\x80\xc1\xd3\ xd3\xd0\xcb\xcc\xcc\xd0\xc0\x86\xbe\xd0\xce\xd1\xce\xd5\xce\xd0\xc4\x88\xc5\xd0\ xce\xc9\xc9\xd4\xd1\xd0\xc3\x85\xd0\xce\xd2\xc9\xcc\xcc\xd2\xd0\xcd\xa8\xcb' >>>
再比较两个库的写入，写入数据gdal用的是WriteRaster。同样DataSet一个，Band一个。另外Band买一送一，还有个WriteArray可用
1 >>> help(ds.WriteRaster)
2 Help on method WriteRaster in module gdal:
3
4 WriteRaster(self, xoff, yoff, xsize, ysize, buf_string, buf_xsize=None, buf_ysiz
5 e=None, buf_type=None, band_list=None) method of gdal.Dataset instance
6
7 >>> help(band.WriteRaster)
8 Help on method WriteRaster in module gdal:
9
10 WriteRaster(self, xoff, yoff, xsize, ysize, buf_string, buf_xsize=None, buf_ysiz
11 e=None, buf_type=None) method of gdal.Band instance
12
13 >>> help(band.WriteArray)
14 Help on method WriteArray in module gdal:
15
16 WriteArray(self, array, xoff=0, yoff=0) method of gdal.Band instance
17
18 >>>
这个就不用做试验了吧。就是把一个二进制字符串放到buf_string位置，那个xoff,yoff,xsize,ysize就是要把数据贴到整张图的什么位置，buf_xsize,buf_ysize是那个二进制字符串表示的区域大小（建议和xsize,ysize一致，不然gdal会自动缩放，而且自动缩放的方法是我们不能控制的，只会用最临近法）。buf_type说明的是输入的二进制字符串的数据类型，如果和原底图的数据类型不一样，会自动扩充和截断数据长度，band_list就是要输入的波段的列表。
PIL中对数据的写入用的是paste，
手册里的解释更多。但差不多的解释就是把im的数据贴到所属对象的box框里。
好了。既然上面验证了两个库中读取的数据二进制一样，那么就可以互相交换数据了。可以把gdal的数据读出，贴到pil打开的数据中，也可以把PIL的数据读出，贴到gdal打开的数据中，尽情地干点偷梁换柱，移花接木的勾当（这个我拿手:-)）。
在PIL中还有一个很好的东东－－fromstring
更多的要去看手册14和16页。Image本身有一个，im（Image打开的对象）也有一个。好用啊！！！创建一个空数据，然后用个fromstring，然后保存，就是一张新图！保存一下刚才那个10*10的图，看看是什么……
另一种:
哈，只要一行！
我喜欢PIL的保存！真是太方便啦！
这么简单的方法，要把GeoTiff转出其他常见格式（不考虑空间信息的时候），就不要用gdal的创建方法了(真是太麻烦了！)用PIL的fromstring再save一下就好了，也不用创建驱动，然后区分Creat和CreateCopy的方式等等步骤，gdal就这点不好……不过如果要转出遥感特殊格式，比如hdf4啊什么的，或者转出的波段比RGB多，还是用gdal的好了。
当然，这里尤其需要注意两个库默认读写哲学的不同。必要的时候要数组转换，或者打开PIL的编码方式。不然转出的东西就不是个东西了。
反馈
如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆这个论坛。