这篇博客源于一个实际需求，在之前利用爬虫方式将Google卫星地图爬下来的时候，虽然可以根据WMTS规则算出每个像素所对应的经纬度，但这个地理信息并不是嵌在图像里的，这样并不利于之后的应用。 因此需要将影像内容和地理信息以GeoTIFF的形式保存下来方便使用。这便是本篇博客主要探讨的问题。之前写过的相关博客列举如下，方便查阅：
之前的这篇博客中介绍了Python中GDAL库的安装与基本使用。 在这篇博客中则介绍了GDAL读写遥感影像的基本方法，本篇博客部分代码以此为基础。 关于WMTS瓦片的一些知识在这篇博客中有介绍。 最后在这篇以及这篇博客中具体介绍了瓦片下载的具体方法。
1.GeoTIFF简介
关于TIFF文件网上的介绍有很多，对于普通应用而言，就知道它是一种图像格式就行了，和jpg、png等等一样。在遥感领域基本所有的影像都是以tif的格式存储的。
但因为专业属性，除了影像内容本身，还会存储一些地理位置、投影信息等信息。
带有这些地理信息的TIFF文件即可以认为是GeoTIFF。
它是一种文件标准，官网是这个。官网中也给出了开源的libgeotiff库用于对其进行读写等基本操作。
但在这篇博客中我们其实不需要过多了解GeoTIFF到底是怎么存储地理信息的，因为我们会使用GDAL库直接生成GeoTIFF文件。所以我们只需要关心GDAL库中如何将地理信息保存到文件中。
具体说来在GDAL库中有GetGeoTransform()与SetGeoTransform()两个函数分别用来读取以及写入地理信息。对于Python的API而言，其输入和输出的是一个包含6个数字的元组，表达一个地理变换。
对于一个geoTransform元组，其各元素表达含义如下：
- geoTransform[0]：左上角像素经度
- geoTransform[1]：影像宽度方向上的分辨率(经度范围/像素个数)
- geoTransform[2]：旋转, 0表示上面为北方
- geoTransform[3]：左上角像素纬度
- geoTransform[4]：旋转, 0表示上面为北方
- geoTransform[5]：影像宽度方向上的分辨率(纬度范围/像素个数)
此外还有SetProjection()和GetProjection()函数用于读取和设置影像的投影信息(坐标系)。投影信息在GDAL中其实就是一串固定格式(如OGC-WKT格式)的字符串，具体可以参见这篇、这篇以及这篇博客，介绍的都很详细。
例如地理坐标(经纬度)常用的WGS84系的WKT描述如下：
GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG", "9122"]], AUTHORITY["EPSG", "4326"]]
具体各项含义见上述博客。在实际使用时，将上述字符串传入即可。 如果不设置坐标系的话，程序也不会报错，但其坐标数值就可能没有意义了(例如无法知道是哪个参考系下的东经120度了)。当然如果真的不需要设置不会产生歧义，直接不用管或者调用这个函数，但传入空字符串就OK了。
2.代码实现与实例
在上面简单介绍后可知，要想将地理信息存入影像，其实只需要在保存影像时调用SetGeoTransform()函数并传入对应的6个变换数字即可。因此代码如下：
# coding=utf-8
from osgeo import gdal
from gdalconst import *
def readImage(img_path):
data = []
# 以只读方式打开遥感影像
dataset = gdal.Open(img_path, GA_ReadOnly)
if dataset is None:
print("Unable to open image file.")
return data
else:
print("Open image file success.")
geoTransform = dataset.GetGeoTransform()
print geoTransform
im_proj = dataset.GetProjection() # 获取投影信息
print im_proj
bands_num = dataset.RasterCount
print("Image height:" + dataset.RasterYSize.__str__() + " Image width:" + dataset.RasterXSize.__str__())
print(bands_num.__str__() + " bands in total.")
for i in range(bands_num):
# 获取影像的第i+1个波段
band_i = dataset.GetRasterBand(i + 1)
# 读取第i+1个波段数据
band_data = band_i.ReadAsArray(0, 0, band_i.XSize, band_i.YSize)
data.append(band_data)
print("band " + (i + 1).__str__() + " read success.")
return data
def writeImage(bands, path, geotrans=None, proj=None):
projection = [
# WGS84坐标系(EPSG:4326)
"""GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]], UNIT["degree", 0.01745329251994328, AUTHORITY["EPSG", "9122"]], AUTHORITY["EPSG", "4326"]]""",
# Pseudo-Mercator、球形墨卡托或Web墨卡托(EPSG:3857)
"""PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]"""
]
if bands is None or bands.__len__() == 0:
return
else:
# 认为各波段大小相等，所以以第一波段信息作为保存
band1 = bands[0]
# 设置影像保存大小、波段数
img_width = band1.shape[1]
img_height = band1.shape[0]
num_bands = bands.__len__()
# 设置保存影像的数据类型
if 'int8' in band1.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in band1.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, img_width, img_height, num_bands, datatype)
if dataset is not None:
if geotrans is not None:
dataset.SetGeoTransform(geotrans) # 写入仿射变换参数
if proj is not None:
if proj is 'WGS84' or proj is 'wgs84' or proj is 'EPSG:4326' or proj is 'EPSG-4326' or proj is '4326':
dataset.SetProjection(projection[0]) # 写入投影
elif proj is 'EPSG:3857' or proj is 'EPSG-3857' or proj is '3857':
dataset.SetProjection(projection[1]) # 写入投影
else:
dataset.SetProjection(proj) # 写入投影
for i in range(bands.__len__()):
dataset.GetRasterBand(i + 1).WriteArray(bands[i])
print("save image success.")
if __name__ == '__main__':
readImage("img.tif")
如下图所示为安徽省蚌埠市淮河文化广场的Google卫星影像图。
利用上述代码可以获取到其地理信息如下：
利用ENVI打开可以看到，能够正确识别出坐标系以及对应的经纬度。
项目及数据放在Github上，点击查看。
3.更为复杂的实例
在一开始提到了，利用爬虫获得的卫星数据是没有地理信息的，因此基于之前已经写好的瓦片爬虫代码以及本篇博客的地理信息相关内容，即可以实现二者的融合，最终生成带有地理信息的、可用的遥感影像。 具体包含三步：下载瓦片、拼接大图、计算并写入地理信息。相关代码及项目传到了Github，点击查看。经测试可以满足预期需求。
4.参考资料
- [1]https://blog.csdn.net/lijie45655/article/details/6771524
本文作者原创,未经许可不得转载,谢谢配合
