博文
“转载”GDAL库学习笔记(七): GDAL和PIL的互操作
||
1. gdal和PIL的互操作
前面提到过，gdal和PIL很相像。它们处理和操作的对象都是栅格图像。但它们又不一样。gdal主要重点放在地理或遥感数据的读写和数据建模以及地理定位和转换，但是PIL的重点是放在图像本身处理上的。
至于在底层数据处理上，两者都可以用Numeric转化的二进制作为数据处理。所以，理论上是可以互相共享和交换数据的。实际上也确实可以。
首先，我要说明的是gdal的核心在波段（band），一切操作的基础和核心都在波段。波段可以单独拿出来操作，至于波段在数据集中的顺序无关紧要。因为遥感图像大多比RGB图像的波段要多，而每个波段单独都是一个完整的整体，每个波段单独拿出来就是一个数据集。而PIL的核心在数据集（DataSet）这里数据集的概念是对应gdal中的数据集的概念，可能在PIL本身中没有这种说法。也就是不把波段单独操作。操作大部分需要RGB一体化地进行，虽然可以单波段操作，但是我想一般人很少这样干！
两部分的操作的主要衔接部分就是创建和读取，以及写入。至于进去了里面怎么变化就是两个库各自的事情了。所以这篇文章的主要内容就是介绍两个库各自的创建，读取和写入的操作，以及两个库的过渡。
比较两个库的读取，gdal读取一个图像中的数据
Toggle line numbers Toggle line numbers
1 >>> import gdal
2 >>> dataset = gdal.Open("E:/gisdata/gtif/sd.tif")
3 >>> ds = dataset
打开了数据集，就有两种方法来获取数据。一种是通过数据集
Toggle line numbers Toggle line numbers
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()
'xd4xd8xd4xcfxccxcexcexc3x7fxbfxd1xd5xcaxd4xd3xd3xd1xc0|xc5x
cdxcdxd1xd1xcaxd0xcexbfx88xc3xd0xd0xcexcbxcaxcaxccxb8wxb9xcc
xcfxcaxcaxd4xd4xcbxb5x87xbexd3xd4xd0xcdxcexbfxcfxc6x80xc1xd3
xd3xd0xcbxccxccxd0xc0x86xbexd0xcexd1xcexd5xcexd0xc4x88xc5xd0
xcexc9xc9xd4xd1xd0xc3x85xd0xcexd2xc9xccxccxd2xd0xcdxa8xcbxce
xd2xcfxcaxc7xcbxcexc5x83xc1xcbxcfxc7xd1xd0xd1xd1xc2x80xc7xc7
xc7xcfxcfxc8xcexcexc1x8cxc5xcaxcaxccxc8xc8xc8xccxba{xbbxc8xc9
xc7xc2xd0xd2xc9xb5x88xc0xcfxcexcdxc5xcaxbdxccxc6x82xc3xd0xcd
xcdxc3xc8xc9xcdxc0x88xc0xcaxc8xcbxc8xcfxcaxccxc0x84xc1xcaxc8
xc5xc5xd0xcdxccxbfx81xccxc8xccxc7xcaxcaxcfxcdxc9xa4xc8xd3xd7
xcdxc8xc5xcdxd3xcdx8dxc9xd0xd4xc7xd1xd0xd5xd6xcax8axcfxccxcc
xd1xd1xcaxd2xd3xc9x96xcdxcfxcfxcexcaxcaxccxd1xc2x85xc3xcdxce
xc9xc8xd6xd6xcexbex93xc6xd4xd3xcfxcbxd0xc2xd4xcdx8axc9xd2xd0
xcfxc9xcexd3xd5xc7x8cxc6xcbxc9xd0xcdxd4xcfxd2xc8x8cxc7xcbxc9
xcaxcaxd5xd2xd1xc5x89xd1xc9xcdxcbxcexcexd1xcfxcexaaxca'
>>> ds.ReadRaster(130,90,10,10)
'xd4xd8xd4xcfxccxcexcexc3x7fxbfxd1xd5xcaxd4xd3xd3xd1xc0|xc5x
cdxcdxd1xd1xcaxd0xcexbfx88xc3xd0xd0xcexcbxcaxcaxccxb8wxb9xcc
xcfxcaxcaxd4xd4xcbxb5x87xbexd3xd4xd0xcdxcexbfxcfxc6x80xc1xd3
xd3xd0xcbxccxccxd0xc0x86xbexd0xcexd1xcexd5xcexd0xc4x88xc5xd0
xcexc9xc9xd4xd1xd0xc3x85xd0xcexd2xc9xccxccxd2xd0xcdxa8xcbxce
xd2xcfxcaxc7xcbxcexc5x83xc1xcbxcfxc7xd1xd0xd1xd1xc2x80xc7xc7
xc7xcfxcfxc8xcexcexc1x8cxc5xcaxcaxccxc8xc8xc8xccxba{xbbxc8xc9
xc7xc2xd0xd2xc9xb5x88xc0xcfxcexcdxc5xcaxbdxccxc6x82xc3xd0xcd
xcdxc3xc8xc9xcdxc0x88xc0xcaxc8xcbxc8xcfxcaxccxc0x84xc1xcaxc8
xc5xc5xd0xcdxccxbfx81xccxc8xccxc7xcaxcaxcfxcdxc9xa4xc8xd3xd7
xcdxc8xc5xcdxd3xcdx8dxc9xd0xd4xc7xd1xd0xd5xd6xcax8axcfxccxcc
xd1xd1xcaxd2xd3xc9x96xcdxcfxcfxcexcaxcaxccxd1xc2x85xc3xcdxce
xc9xc8xd6xd6xcexbex93xc6xd4xd3xcfxcbxd0xc2xd4xcdx8axc9xd2xd0
xcfxc9xcexd3xd5xc7x8cxc6xcbxc9xd0xcdxd4xcfxd2xc8x8cxc7xcbxc9
xcaxcaxd5xd2xd1xc5x89xd1xc9xcdxcbxcexcexd1xcfxcexaaxca'
>>>
--------------------------------------------------------------------------------
从波段中获取数据和从数据集中获取数据有十分相似的方法。
Toggle line numbers Toggle line numbers
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打开获取数据的过程如下：
Toggle line numbers Toggle line numbers
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()
'xd4xcexd3xd8xd2xd7xd4xcfxcdxcfxcaxc8xccxc7xc5xcexcbxcdxcexc
exd3xc3xc5xcdx7fx83x8dxbfxc1xc9xd1xcbxd0xd5xcfxd4xcaxc7xc7xd
4xd1xd1xd3xd0xd0xd3xd1xd5xd1xd1xd6xc0xc2xca|x80x8axc5xc7xcfx
cdxc7xccxcdxc7xccxd1xcfxd1xd1xcfxd1xcaxc8xcaxd0xcexd2xcexcex
d3xbfxc1xc9x88x8cx96xc3xc5xcdxd0xcaxcfxd0xcaxcfxcexccxcexcbx
c8xcaxcaxc8xcaxcaxc8xccxccxccxd1xb8xbaxc2w{x85xb9xbbxc3xccxc8
xcdxcfxc9xcexcaxc7xc9xcaxc2xc8xd4xd0xd6xd4xd2xd6xcbxc9xcexb5
xb5xbex87x88x93xbexc0xc6xd3xcfxd4xd4xcexd3xd0xcdxcfxcdxc5xcb
xcexcaxd0xbfxbdxc2xcfxccxd4xc6xc6xcdx80x82x8axc1xc3xc9xd3xd0
xd2xd3xcdxd0xd0xcdxcfxcbxc3xc9xccxc8xcexccxc9xd3xd0xcdxd5xc0
xc0xc7x86x88x8cxbexc0xc6xd0xcaxcbxcexc8xc9xd1xcbxd0xcexc8xcd
xd5xcfxd4xcexcaxcfxd0xccxd2xc4xc0xc8x88x84x8cxc5xc1xc7xd0xca
xcbxcexc8xc9xc9xc5xcaxc9xc5xcaxd4xd0xd5xd1xcdxd2xd0xccxd1xc3
xbfxc5x85x81x89xd0xccxd1xcexc8xc9xd2xccxcdxc9xc7xcbxccxcaxce
xccxcaxcexd2xcfxd1xd0xcdxcfxcdxc9xcexa8xa4xaaxcbxc8xca'
>>>
这里注意，在pil的下，截取区域矩形的定义和gdal不同，gdal是顶点X、顶点Y、宽、高；pil是顶点X、顶点Y、终点X，终点Y
诶，似乎不一样啊！
确实不一样，不过这就是gdal和pil的区别了。转换一下吧。
Toggle line numbers Toggle line numbers
1 >>> data = ds.ReadAsArray(130,90,10,10)
2 >>> datas = [i for i in data]
3 >>> from Numeric import *
4 >>> datas = [reshape(i,(-1,1)) for i in data]
5 >>> datas = concatenate(datas,1)
6 >>>
激动人心的时刻到了！当当当当～～～～～
--------------------------------------------------------------------------------
>>> datas.tostring()
'xd4xcexd3xd8xd2xd7xd4xcfxcdxcfxcaxc8xccxc7xc5xcexcbxcdxcexc
exd3xc3xc5xcdx7fx83x8dxbfxc1xc9xd1xcbxd0xd5xcfxd4xcaxc7xc7xd
4xd1xd1xd3xd0xd0xd3xd1xd5xd1xd1xd6xc0xc2xca|x80x8axc5xc7xcfx
cdxc7xccxcdxc7xccxd1xcfxd1xd1xcfxd1xcaxc8xcaxd0xcexd2xcexcex
d3xbfxc1xc9x88x8cx96xc3xc5xcdxd0xcaxcfxd0xcaxcfxcexccxcexcbx
c8xcaxcaxc8xcaxcaxc8xccxccxccxd1xb8xbaxc2w{x85xb9xbbxc3xccxc8
xcdxcfxc9xcexcaxc7xc9xcaxc2xc8xd4xd0xd6xd4xd2xd6xcbxc9xcexb5
xb5xbex87x88x93xbexc0xc6xd3xcfxd4xd4xcexd3xd0xcdxcfxcdxc5xcb
xcexcaxd0xbfxbdxc2xcfxccxd4xc6xc6xcdx80x82x8axc1xc3xc9xd3xd0
xd2xd3xcdxd0xd0xcdxcfxcbxc3xc9xccxc8xcexccxc9xd3xd0xcdxd5xc0
xc0xc7x86x88x8cxbexc0xc6xd0xcaxcbxcexc8xc9xd1xcbxd0xcexc8xcd
xd5xcfxd4xcexcaxcfxd0xccxd2xc4xc0xc8x88x84x8cxc5xc1xc7xd0xca
xcbxcexc8xc9xc9xc5xcaxc9xc5xcaxd4xd0xd5xd1xcdxd2xd0xccxd1xc3
xbfxc5x85x81x89xd0xccxd1xcexc8xc9xd2xccxcdxc9xc7xcbxccxcaxce
xccxcaxcexd2xcfxd1xd0xcdxcfxcdxc9xcexa8xa4xaaxcbxc8xca'
>>>
--------------------------------------------------------------------------------
相同了吧！
这里就表现了两个库的设计哲学不同了。gdal的核心是band，读取的数据是默认以band组织的，pil的核心是dataset，读取的数据默认是统一组织的。换句话说，gdal看{{RRR……}{GGG……}{BBB……}}比较顺眼，pil看{RGB,RGB,RGB……}比较顺眼。
同时把数据降低到Band水平（其实不能说降低，Band的水平不比DataSet低），就可以看出来在Band水平，两个库是一样一样一样的啊！
--------------------------------------------------------------------------------
>>> r,g,b = region.split()
>>> r.tostring()
'xd4xd8xd4xcfxccxcexcexc3x7fxbfxd1xd5xcaxd4xd3xd3xd1xc0|xc5x
cdxcdxd1xd1xcaxd0xcexbfx88xc3xd0xd0xcexcbxcaxcaxccxb8wxb9xcc
xcfxcaxcaxd4xd4xcbxb5x87xbexd3xd4xd0xcdxcexbfxcfxc6x80xc1xd3
xd3xd0xcbxccxccxd0xc0x86xbexd0xcexd1xcexd5xcexd0xc4x88xc5xd0
xcexc9xc9xd4xd1xd0xc3x85xd0xcexd2xc9xccxccxd2xd0xcdxa8xcb'
>>> band = ds.GetRasterBand(1)
>>> band.ReadRaster(130,90,10,10)
'xd4xd8xd4xcfxccxcexcexc3x7fxbfxd1xd5xcaxd4xd3xd3xd1xc0|xc5x
cdxcdxd1xd1xcaxd0xcexbfx88xc3xd0xd0xcexcbxcaxcaxccxb8wxb9xcc
xcfxcaxcaxd4xd4xcbxb5x87xbexd3xd4xd0xcdxcexbfxcfxc6x80xc1xd3
xd3xd0xcbxccxccxd0xc0x86xbexd0xcexd1xcexd5xcexd0xc4x88xc5xd0
xcexc9xc9xd4xd1xd0xc3x85xd0xcexd2xc9xccxccxd2xd0xcdxa8xcb'
>>>
--------------------------------------------------------------------------------
再比较两个库的写入，写入数据gdal用的是WriteRaster。同样DataSet一个，Band一个。另外Band买一送一，还有个WriteArray可用
Toggle line numbers Toggle line numbers
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，
Toggle line numbers Toggle line numbers
1 >>> help(im.paste)
2 Help on method paste in module Image:
3
4 paste(self, im, box=None, mask=None) method of TiffImagePlugin.TiffImageFile ins
5 tance
6 Paste other image into region
7
8 >>>
手册里的解释更多。但差不多的解释就是把im的数据贴到所属对象的box框里。
好了。既然上面验证了两个库中读取的数据二进制一样，那么就可以互相交换数据了。可以把gdal的数据读出，贴到pil打开的数据中，也可以把PIL的数据读出，贴到gdal打开的数据中，尽情地干点偷梁换柱，移花接木的勾当（这个我拿手:-)）。
在PIL中还有一个很好的东东－－fromstring
Toggle line numbers Toggle line numbers
1 >>> help(Image.fromstring)
2 Help on function fromstring in module Image:
3
4 fromstring(mode, size, data, decoder_name='raw', *args)
5 Load image from string
6
7 >>>
更多的要去看手册14和16页。Image本身有一个，im（Image打开的对象）也有一个。好用啊！！！创建一个空数据，然后用个fromstring，然后保存，就是一张新图！保存一下刚才那个10*10的图，看看是什么……
Toggle line numbers Toggle line numbers
1 >>> newim = Image.new("RGB",(10,10))
2 >>> newim.fromstring(datas.tostring())
3 >>> newim.save("f:/png/new.png","png")
4 >>>
另一种:
Toggle line numbers Toggle line numbers
1 >>> Image.fromstring("RGB",(10,10),datas.tostring()).save("f:/png/new1.png","p
2 ng")
3 >>>
哈，只要一行！
我喜欢PIL的保存！真是太方便啦！
这么简单的方法，要把GeoTiff转出其他常见格式（不考虑空间信息的时候），就不要用gdal的创建方法了(真是太麻烦了！)用PIL的fromstring再save一下就好了，也不用创建驱动，然后区分Creat和CreateCopy的方式等等步骤，gdal就这点不好……不过如果要转出遥感特殊格式，比如hdf4啊什么的，或者转出的波段比RGB多，还是用gdal的好了。
当然，这里尤其需要注意两个库默认读写哲学的不同。必要的时候要数组转换，或者打开PIL的编码方式。不然转出的东西就不是个东西了。
本篇文章来源于 GIS空间站 转载请以链接形式注明出处 网址：http://www.gissky.net/Article/649_3.htm
https://blog.sciencenet.cn/blog-106723-221145.html
上一篇：与测量相关的资源网站
下一篇：TGO使用说明书
