我要投稿
投稿指南
RSS订阅
3. 我的理解
上面文档说得很玄，还是让我们看看坐标系统究竟是什么样子。
1 >>> import gdal
2 >>> dataset = gdal.Open("e:/gisdata/gtif/spot.tif")
3 >>> dir(dataset)
4 ['AddBand', 'AdviseRead', 'BuildOverviews', 'FlushCache', 'GetDescription', 'Get
5 Driver', 'GetGCPCount', 'GetGCPProjection', 'GetGCPs', 'GetGeoTransform', 'GetMe
6 tadata', 'GetProjection', 'GetProjectionRef', 'GetRasterBand', 'GetSubDatasets',
7 'RasterCount', 'RasterXSize', 'RasterYSize', 'ReadAsArray', 'ReadRaster', 'Refr
8 eshBandInfo', 'SetDescription', 'SetGCPs', 'SetGeoTransform', 'SetMetadata', 'Se
9 tProjection', 'WriteRaster', '__del__', '__doc__', '__init__', '__module__', '_b
10 and', '_o']
11 >>> dataset.GetProjectionRef()
12 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
13 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
14 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
15 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
16 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
17 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
18 SG","9001"]],AUTHORITY["EPSG","26713"]]'
19 >>>
ERROR: EOF in multi-line statement
3.1. WKT
最后那一大串字符串就代表了坐标系统了，没有缩进编排会看晕了（这个东西我第一次看以为是脚本）。GDAL的坐标系统的标示方法是用OpenGIS的WKT字符串表示。当然也混合了EPSG权威编码，关于WKT可以看这里 。 它是绝对的标准，完全解释了坐标系统的组成，结构和层次，它是由一些层次嵌套构成的。
1 Coordinate System WKT
2 <coordinate system> = <horz cs> | <geocentric cs> | <vert cs> | <compd cs> | <fitted cs> | <local cs>
3 <horz cs> = <geographic cs> | <projected cs>
4 <projected cs> = PROJCS["" , <geographic cs>, <projection>, {<parameter>,}* <linear unit> {,<twin axes>}{,<authority>}]
5 <projection> = PROJECTION["" {,<authority>}]
6 <geographic cs> = GEOGCS["" , <datum>, <prime meridian>, <angular unit> {,<twin axes>} {,<authority>}]
7 <datum> = DATUM["" , <spheroid> {,<to wgs84>} {,<authority>}]
8 <spheroid> = SPHEROID["" , <semi-major axis>, <inverse flattening> {,<authority>}]
9 <semi-major axis> = <number>
10 <inverse flattening> = <number>
11 <prime meridian> = PRIMEM["" , <longitude> {,<authority>}]
12 <longitude> = <number>
13 <angular unit> = <unit>
14 <linear unit> = <unit>
15 <unit> = UNIT["" , <conversion factor> {,<authority>}]
16 <conversion factor> = <number>
17 <geocentric cs> = GEOCCS["" , <datum>, <prime meridian>, <linear unit> {,<axis>, <axis>, <axis>} {,<authority>}]
18 <authority> = AUTHORITY["" , ""]
19 <vert cs> = VERT_CS["" , <vert datum>, <linear unit>, {<axis>,} {,<authority>}]
20 <vert datum> = VERT_DATUM["" , <datum type> {,<authority>}]
21 <datum type> = <number>
22 <compd cs> = COMPD_CS["" , <head cs>, <tail cs> {,<authority>}]
23 <head cs> = <coordinate system>
24 <tail cs> = <coordinate system>
25 <twin axes> = <axis>, <axis>
26 <axis> = AXIS["" , NORTH | SOUTH | EAST | WEST | UP | DOWN | OTHER]
27 <to wgs84s> = TOWGS84[<seven param>]
28 <seven param> = <dx>, <dy>, <dz>, <ex>, <ey>, <ez>, <ppm>
29 <dx> = <number>
30 <dy> = <number>
31 <dz> = <number>
32 <ex> = <number>
33 <ey> = <number>
34 <ez> = <number>
35 <ppm> = <number>
36 <fitted cs> = FITTED_CS["" , <to base>, <base cs>]
37 <to base> = <math transform>
38 <base cs> = <coordinate system>
39 <local cs> = LOCAL_CS["" , <local datum>, <unit>, <axis>, {,<axis>}* {,<authority>}]
40 <local datum> = LOCAL_DATUM["" , <datum type> {,<authority>}]
这个结构的阅读方法是这样的，所有的尖括号<>中的值都是一个需要被替代的东西（我觉得应该叫做一个类型）。你可以在表的其余中找到对应的类型描述，但是这个类型是不能直接写出来的，也就是说相当于一个链接，也可以看成是一个指针或一个引用，它指向一个对象，但是却不是一个对象，你需要用一个确实的东西来代替它。而且文档的下面有各种类型的解释。在花括号{}中的是可选的项目，表现在实际的坐标表示中，就是可写可不写。等号=是用来表示类型的表达形式。 用 | 分割开的是并列的选项，也就是在实际坐标表示中你只能选用其中的一个。*号表示有数个，多个相同的类型组成。除了上面的其他的符号就代表确实的东西（[]是需要实际写入坐标系表示字符串的一部分，不是这个语法表达表中的表达符号，冒号" 逗号, 同样也是实体的一部分），是需要到坐标系中的。最后注意书写顺序从 解释起来很麻烦(汗都说出来了)，如果你还是不理解，算了，就当我没说，其实你可以完全不管这些的（但是如果你要写，可以看看这里，这里来获取更多有关于坐标系统和其参数的信息，你还可以下载postgresql，如果安装了postgis插件，你就拥有了几乎所有常用的坐标系统的wkt表示）。如果我们不亲自写不同坐标间的转换函数。我们只要关心两个坐标系统是否完全一样，其他都可以不去管它。要是说起不同坐标系统之间转换，又有一大堆的东西要说了，还是以后开个专门的页面写写吧。对于我们来说重要的是另外一个部分－－坐标转换。
这是我们主要要对付的。只有通过这几个数据，我们才可以知道自己手中这张图该如何铺展到绘制平面上。上面说过了，栅格坐标到地理坐标之间的转换可以通过两种方式来进行：一种是仿射坐标转换，一种是GCPs。一般来说，仿射坐标转换比较常见。仿射坐标转换就是通过几个值来进行栅格框架到地理框架的映射。上面的例子中出现了6个值。这六个值可以用下面的
公式来解释，而且数组顺序对应起来了，理解起来也不困难。至于几个参数究竟是什么意思，这里讲得更清楚：
看起来和ESRI为GeoTiff定义的tfw文件格式很像。至于是否真的是完全对应，我还不敢说。
有了上面的公式，我们就可以把栅格的每个点代入公式来求出其在地球上的实际位置（当然，求出的数值的单位是在坐标系统中定义的那个，而非经纬度坐标，要求经纬度坐标，还要通过一步转换）。当然我们可以不用这么麻烦，只要求出整个图像外框上下左右四个边界就可以定下整个图像的位置，因为图像都是矩形的（那些非矩形的图像其实是有另外一层二值的掩膜图像（mask）叠加而成的，这种图像只有两种值，一种是0,一种是1，1代表按照图像原色显示，0代表不显示）。图像中某点的地理位置可以通过外框四点线性计算。
现在看看我们这张图的边界地理范围：
上边框和左边框都不用算了，是4928000和590000。 求右下角的坐标，代入公式：
右边框 = Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2) = 590000.0 + 950*20.0 + 700*0.0 = 609000.0
下边框 = Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5) = 4928000.0 + 950*0.0 + 700*(-20.0) = 4914000.0
开ArcGIS验证一下，没错。
好了，有了这些，就可以把图像和地理坐标联系起来了。有了地理坐标，就可以把图像在照屏幕上预先定好的坐标系画出来了。
需要注意的一点是，对于GeoTiff数据源来说，坐标系统和仿射变换参数有两种存储方式，一种是直接存储在.tif文件内部，这种存储方式是按照tif内部的的键值对方式来存储的。具体的一些细节可以参看这里。官方的说明pdf可以在这里下载。也可以把图像仅仅存储到tif文件中，而把坐标系统和仿射参数什么的单独提取出来，创建两个文件，一个是prj文件，存放WKT坐标系字符串，一个是tfw文件，存放仿射转换参数。
是对应的实际的类型的名称和编码。不是在表中有对应。就如上面的例子，基准面名是"North_American_Datum_1927"，对应的3.2. 坐标转换
1 >>> dataset.GetGCPs()
2 []
3 >>> dataset.GetGeoTransform()
4 (590000.0, 20.0, 0.0, 4928000.0, 0.0, -20.0)
5 >>>
1 Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
2 Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
1 adfGeoTransform[0] /* top left x */
2 adfGeoTransform[1] /* w-e pixel resolution */
3 adfGeoTransform[2] /* rotation, 0 if image is "north up" */
4 adfGeoTransform[3] /* top left y */
5 adfGeoTransform[4] /* rotation, 0 if image is "north up" */
6 adfGeoTransform[5] /* n-s pixel resolution */
1 >>> dataset.GetGeoTransform()
2 (590000.0, 20.0, 0.0, 4928000.0, 0.0, -20.0)
3 >>> dataset.RasterXSize
4 950
5 >>> dataset.RasterYSize
6 700
7 >>>
1 In projected or local coordinates
2 Left: 590000.000000
3 Right: 609000.000000
4 Top: 4928000.000000
5 Bottom: 4914000.000000
