前面提到过,GDAL和PIL很相像。它们处理和操作的对象都是栅格图像。 但它们又不一样。 GDAL主要重点放在地理或遥感数据的读写和数据建模以及地理定位和转换, 但是PIL的重点是放在图像本身处理上的。
至于在底层数据处理上,两者都可以用 numpy 转化的二进制作为数据处理。 所以,理论上是可以互相共享和交换数据的。实际上也确实可以。
首先,我要说明的是GDAL的核心在波段(band), 一切操作的基础和核心都在波段。 波段可以单独拿出来操作,至于波段在数据集中的顺序无关紧要。 因为遥感图像大多比RGB图像的波段要多,而每个波段单独都是一个完整的整体, 每个波段单独拿出来都是一个数据集。而 Pillow 的核心在数据集(DataSet)这里的概念是对应GDAL中的数据集的概念。 当然,在 Pillow 本身中没有这种说法,也就是不把波段单独操作, 操作大部分需要RGB一体化地进行。
两部分的操作的主要衔接部分就是创建、读取与写入。 读取数据后怎么处理是两个库各自的事情。 所以这里主要内容就是介绍两个库各自的创建,读取和写入的操作,以及两个库的过渡。
from osgeo import gdal
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
data_arr = dataset.ReadAsArray(30,70,5,5)
type(data_arr)
/opt/conda/lib/python3.12/site-packages/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default. warnings.warn(
numpy.ndarray
data_arr
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)
data_bin = dataset.ReadRaster(30,70,5,5)
data_bin
bytearray(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读出的相同。
data_arr.tobytes()
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'
从波段中获取数据和从数据集中获取数据的方法十分相似。
from PIL import Image
im = Image.open("/data/gdata/geotiff_file.tif")
region = im.crop((30,70,35,75))
region.tobytes()
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()
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.tobytes()
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'
r,g,b = region.split()
r.tobytes()
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'
band = dataset.GetRasterBand(1)
band.ReadRaster(30,70,5,5)
bytearray(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')