Tiff与gdal

21 Nov 2018

1.tiff格式

标签图像文件格式（Tag Image File Format，简写为TIFF）是一种灵活的位图格式，主要用来存储包括照片和艺术图在内的图像。它最初由Aldus公司与微软公司一起为PostScript打印开发。TIFF与JPEG和PNG一起成为流行的高位彩色图像格式。

TIFF 是一个灵活适应性强的文件格式，通过在文件头中包含“标签”它能够在一个文件中处理多幅图像和数据。标签能够标明图像的如图像大小这样的基本几何尺寸或者定义图像数据是如何排列的并且是否使用了各种各样的图像压缩选项。

TIFF有一个使用LZW压缩的选项，这是一种减小文件大小的无损技术。与JPEG不同，TIFF文件可以编辑然后重新存储而不会有压缩损失。

TIFF文件以.tif为扩展名。其数据格式是一种3级体系结构，从高到低依次为：文件头、一个或多个称为IFD的包含标记指针的目录和数据。

2.GTiff – GeoTIFF File Format

GeoTIFF利用了TIFF的可扩展性，在其基础上加了一系列标志地理信息的标签（Tag），来描述卫星成像系统、航空摄影、地图信息和DEM等。 GeoTIFF设计使得标准的地图坐标系定义可以随意存储为单一的注册标签。GeoTIFF也支持非标准坐标系的描述，为了在不同的坐标系间转换，可以通过使用3~4个另设的TIFF标签来实现。GeoTIFF 使用了 GeoKeys 来组织众多空间参数，所有 GeoKeys 由 GeoKeyDirectoryTag 来索引。

GeoTIFF 支持三种坐标空间 : 栅格空间 (Raster Space) 、设备空间 (Device Space) 和模型空间 (Model Space) 。栅格空间是存储图像的行列号的坐标系统。它有 2 种形式 : 一是 PixelIs-Area 型, 主要用于非 DEM 的数据 ; 二是 PixelIsPoint 型, 主要用于 DEM 。设备空间是使用 TIFF 格式中定义的 6 个基本 Tag 来描述图像的分辨率单位及图像定位。模型空间是 GeoTIFF 图像的栅格坐标所对应的实际地图的经纬度坐标或直角坐标。

各个地理标签的具体含义如下 :

( 1)ModelPixelScaleTag( 像元比例 ) : 存放着图像中的某一点在栅格空间中的坐标与在模型空间中的坐标的缩放比例。

( 2)ModelTiepointTag( 控制点 ) : 图像中栅格坐标与其对应的模型坐标形成的坐标控制点对。

( 3)ModelTransformationTag( 变换矩阵 ) : 含有 16 个双精度 ( 4)GeoDoubleParamsTag( 双精度参数 ) : 代表 GeoTIFF 定义的一种数据类型 , 用来存储双精度类型的地理键 (GeoKeys) 。

( 5)GeoAsciiParamsTag(ASCII 参数 ) : 存储字符型的地理键 (GeoKey) 值 , 保存字符型的地理键 (GeoKeys) 。

( 6)GeoKeyDirectoryTag( 地理信息目录 ) : 是 6 个地理标签中最重要、最复杂的一个 , 可分为头和记录两部分 : 头部的结构为 Header={ 目录版本号 , 修订版本号 , 副版本号 , 地理键的个数 }; 每条记录的结构为 KeyEntry={ 地理键 ID, 存放位置 , 元素的个数 , 值 / 索引 } 。 ID 号唯一标识了地理键 , 存放位置表示地理键存放在哪个标签中 ( 主要指 GeoDoubleParamsTag 和 GeoAscii-ParamsTag) 。如果值为 0 表示该键为短整型、个数为 1 则它的值就保存在记录中。否则 , 其类型由 TIFFTagLocation 暗指 , 值即存放在 TIFFTagLocation 指定的标签中 , 第一个元素在标签中的索引为偏移量。

When built with internal libtiff or with libtiff >= 4.0, GDAL also supports reading and writing BigTIFF files (evolution of the TIFF format to support files larger than 4 GB).

3.GDAL读取TIFF高程

GDAL是一个操作各种栅格地理数据格式的库。包括读取、写入、转换、处理各种栅格数据格式。它使用了一个单一的抽象数据模型就支持了大多数的栅格数据。除了栅格操作，这个库还同时包括了操作矢量数据的另一个有名的库OGR，这样这个库就同时具备了操作栅格和矢量数据的能力。

GDAL使用抽象数据模型(abstractdatamodel)来解析它所支持的数据格式，抽象数据模型包括数据集(dataset)，坐标系统，仿射地理坐标转换(Affine GeoTransform)，大地控制点(GCPs)，元数据(Metadata)，栅格波段(RasterBand)，颜色表(ColorTable)，子数据集域(Subdatasets Domain)，图像结构域(Image_StructureDomain)，XML域(XML:Domains)。

GDAL基础类

GDALMajorObject类：抽象类，带有元数据的对象。
GDALDdataset类：通常是从一个栅格文件中提取的相关联的栅格波段集合和这些波段的元数据； GDALDdataset也负责所有栅格波段的地理坐标转换(georeferencing transform)和坐标系定义。Dataset的坐标系统由OpenGIS WKT字符串定义。　
GDALDriver类：文件格式驱动类，GDAL会为每一个所支持的文件格式创建一个该类的实体，来管理该文件格式。
GDALDriverManager类：文件格式驱动管理类，用来管理GDALDriver类。

读取某一像素点的值，需要分两步首先读取一个波段(band)：GetRasterBand()，其参数为波段的索引号然后用ReadAsArray(, , , )，读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。如果将矩阵大小设为1X1，就是读取一个像素了。但是这一方法只能将读出的数据放到矩阵中，就算只读取一个像素也是一样。
例如：
band = ds.GetRasterBand(1)
data = band.ReadAsArray(xOffset, yOffset, 1, 1) 如果想一次读取一整张图，那么将offset都设定为0，size则设定为整个图幅的size，

例如：
data = band.ReadAsArray(0, 0, cols, rows) 但是要注意，从data中读取某一像素的值，必须要用data[yoff, xoff]。注意不要搞反了。数学中的矩阵是[row,col]，而这里恰恰相反！这里面row对应y轴，col对应x轴。

利用PYTHON与MATLIB可以查看高程，只是有点卡，毕竟不是专业工具。

dem

仿射地理坐标转换（Affine GeoTransform）

GDAL数据集有两种方式来表示栅格行列号坐标和地理坐标之间的关系，第一种，也是最常用的，就是使用仿射变换来表示（另外一种是GCP点）。

仿射变换包含六个参数，可以使用函数 GDALDataset::GetGeoTransform()获取，图像行列号和地理空间之间的变换关系如下：

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5) 在北方向上的图像中，参数GT(2)和GT(4)的值是0，GT(1)表示象元宽度，GT(5)表示象元高度。点 (GT(0),GT(3))表示图像左上角的横纵坐标值。

值得注意的是：行列号坐标系统中，图像的（0.0，0.0）坐标在图像的左上角，坐标（象元宽个数,象元高个数）表示图像的右下角坐标。图像的左上角象元的中心点坐标是（0.5,0.5）。

4.坐标系

Google Maps、Virtual Earth等网络地理所使用的地图投影，常被称作Web Mercator或Spherical Mercator，它与常规墨卡托投影的主要区别就是把地球模拟为球体而非椭球体。

Web墨卡托投影坐标系：

以整个世界范围，赤道作为标准纬线，本初子午线作为中央经线，两者交点为坐标原点，向东向北为正，向西向南为负。

X轴：由于赤道半径为6378137米，则赤道周长为2PIr = 2*20037508.3427892，因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。

Y轴：由墨卡托投影的公式可知，同时上图也有示意，当纬度φ接近两极，即90°时，y值趋向于无穷。这是那些“懒惰的工程师”就把Y轴的取值范围也限定在[-20037508.3427892,20037508.3427892]之间，搞个正方形。

经度：全球范围：[-180,180]。

纬度：上面已知，纬度不可能到达90°，懒人们为了正方形而取的-20037508.3427892，经过反计算，可得到纬度85.05112877980659。因此纬度取值范围是[-85.05112877980659，85.05112877980659]。其余的地区怎么办？没事，企鹅们不在乎。

因此，地理坐标系（经纬度）对应的范围是：最小(-180,-85.05112877980659)，最大(180, 85.05112877980659)。

#经纬度转Wev墨卡托
def lonLat2WebMercator(lonLat)
    mercator = vec3d()
    x = lonLat.x *20037508.34/180
    y = math.log(math.tan(( 90 + lonLat.y ) * math.pi / 360)) / (math.pi/180)
    y = y *20037508.34/180;
    mercator.x = x;
    mercator.y = y;
    return mercator

#Web墨卡托转经纬度
def WebMercator2lonLat( mercator ):
    lonLat= vec3d()
    x = mercator.x/20037508.34*180
    y = mercator.y/20037508.34*180
    y= 180/ math.pi *( 2* math.atan( math.exp( y * math.pi / 180)) - math.pi / 2 )
    lonLat.x = x
    lonLat.y = y
    return lonLat

取高程数据示例

import gdal
from gdalconst import *
import numpy as np
import math

class vec2d:
    def __init__(self, x=0,y=0):
        self.x=x
        self.y=y

class vec3d:
    def __init__(self, x=0,y=0,z=0):
        self.x=x
        self.y=y
        self.z=z

#经纬度转Wev墨卡托
def lonLat2WebMercator(lonLat)
    mercator = vec3d()
    x = lonLat.x *20037508.34/180
    y = math.log(math.tan(( 90 + lonLat.y ) * math.pi / 360)) / (math.pi/180)
    y = y *20037508.34/180;
    mercator.x = x;
    mercator.y = y;
    return mercator

#Web墨卡托转经纬度
def WebMercator2lonLat( mercator ):
    lonLat= vec3d()
    x = mercator.x/20037508.34*180
    y = mercator.y/20037508.34*180
    y= 180/ math.pi *( 2* math.atan( math.exp( y * math.pi / 180)) - math.pi / 2 )
    lonLat.x = x
    lonLat.y = y
    return lonLat

def getDem():
    gdal.AllRegister()
    dataset = gdal.Open("./L13/sample.tif",GA_ReadOnly)
    if dataset is None:
        print('data is invalid!')
        return -1
    print('driver:', dataset.GetDriver().ShortName)
    print('size: x=',dataset.RasterXSize,
          ',y=',dataset.RasterYSize,',RasterCount=',dataset.RasterCount)#

    im_proj = dataset.GetProjection() #地图投影信息

    nrows = dataset.RasterXSize
    ncols  = dataset.RasterYSize  #这两个行就是读取数据的行列数
    band = dataset.GetRasterBand(3)  #用gdal去读写你的数据，当然dem只有一个波段
    adfGeoTransform = dataset.GetGeoTransform()
    Xmin = adfGeoTransform[0] #你的数据的平面四至
    Ymin = adfGeoTransform[3]
    ori = vec3d()
    ori.x = Xmin; ori.y=Ymin
    ori = WebMercator2lonLat(ori)
    print('LEFT UP X=',ori.x)
    print('LEFT UP Y=',ori.y)
    Xmax = adfGeoTransform[0] + nrows * adfGeoTransform[1] + ncols * adfGeoTransform[2]
    Ymax = adfGeoTransform[3] + nrows * adfGeoTransform[4] + ncols * adfGeoTransform[5]
    x = np.linspace(Xmin,Xmax, ncols)
    y = np.linspace(Ymin,Ymax, nrows)
    z = band.ReadAsArray(0, 0) #这一段就是讲数据的x，y，z化作numpy矩阵

4.DEM高程数据及获取

DEM是数字高程模型的英文简称(Digital Elevation Model)，是研究分析地形、流域、地物识别的重要原始资料。由于DEM 数据能够反映一定分辨率的局部地形特征，因此通过DEM 可提取大量的地表形态信息，可用于绘制等高线、坡度图、坡向图、立体透视图、立体景观图，并应用于制作正射影像、立体地形模型与地图修测。在测绘、水文、气象、地貌、地质、土壤、工程建设、通讯、军事等国民经济和国防建设以及人文和自然科学领域有着广泛的应用。如在工程建设上，可用于如土方量计算、通视分析等；在防洪减灾方面，DEM是进行水文分析如汇水区分析、水系网络分析、降雨分析、蓄洪计算、淹没分析等的基础; 在无线通讯上，可用于蜂窝电话的基站分析等。

介绍下我下载过的数据源：

公开的DEM数据SRTM3 ：地址http://srtm.csi.cgiar.org/srtmdata/,250米精度

https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Africa/，我国境内的数据在Eurasia目录下，每经纬度方格一个文件，文件命名方法是X1X2X3X4.hgt.zip，X1是N或 S表示南北，X2是下方纬度数，X3是E或W表示东西，X4是左方经度数，Documentation目录下有这些数据的介绍

ETOPO

ETOPO1是分辨率为1弧分的全球地形起伏模型，其包含了陆地地形和海洋水深的数据，是目前已知的分辨率最高的地形起伏数据。

其分为两个版本，Ice Surface和Bedrock，两个版本基本一致。不同之处在于在处理南极洲和Greenland地形时，前者给出的是加上冰盖层之后的高程，后者给出的是岩床的高程。

对于每个版本又细分为grid-registered和cell-registered，其中grid-registered是权威版本，cell-registered是衍生版本，因而推荐下载使用grid-registered版本。

https://www.ngdc.noaa.gov/mgg/global/global.html

参考：