GDAL和 Pillow 的互操作

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

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

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

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

使用GDAL读取数据

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

In [2]:
from osgeo import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
data_arr = dataset.ReadAsArray(30,70,5,5)
type(data_arr)
Out[2]:
numpy.ndarray
In [3]:
data_arr
Out[3]:
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 [4]:
data_bin = dataset.ReadRaster(30,70,5,5)
data_bin
Out[4]:
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 [5]:
data_arr.tostring()
Out[5]:
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 [6]:
from PIL import Image
im = Image.open("/gdata/geotiff_file.tif")
In [7]:
region = im.crop((30,70,35,75))
region.tobytes()
Out[7]:
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 [8]:
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[8]:
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 [9]:
r,g,b = region.split()
r.tobytes()  
Out[9]:
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 [10]:
band = dataset.GetRasterBand(1)
band.ReadRaster(30,70,5,5)
Out[10]:
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")

只要一行!


In [11]:
import helper; helper.info()
页面更新时间: 2019-04-14 21:00:33
操作系统/OS: Linux-4.9.0-8-amd64-x86_64-with-debian-9.8
Python: 3.5.3