tfw坐标定义.txt

栅格图像坐标文件TFW中的说明,还是人家的,不知有用没?我觉得是有用的

http://www.gissky.net/Article/qy/200702/646.htm
作者：lilin    文章来源：啄木鸟CPUG站    点击数：4774    更新时间：2007-2-8 
 
1. 关于空间参考
好啦，现在让我们开始进入地理的世界 

一幅地图除了需要画出地物的形状，更重要的是要画出地物的位置。做GIS的头一项就是要做定位。这也就是GIS数据源和普通数据源的区别－－它带有空间信息。GDAL和PIL等图像处理库的最大区别也就在此处。 

空间信息包括两个部分。一个是坐标系统，还有一个就是栅格坐标和地理坐标之间的转换。 

坐标系统是一个很大的概念。包括很多知识点，可以写厚厚的一本书。我们这里也不想多说，毕竟我在这里不是要普及地理知识，我只在这里讲GDAL中涉及的坐标系统概念。因为我自己对坐标系统学习也并不是太深入（我不是测绘专业的），所以，我的一些理解可能是错误的，希望大家能及时指正。 根据GDAL数据模型文档的解释，GDAL数据集类中包括的坐标系统解释是这样的： 


2. GDAL文档的解释
数据集的坐标系统由OpenGIS WKT字符串定义，它包含了： 

1、一个总体的坐标系名 

2、一个地理图形坐标系统名 

3、一个基准面定义 

4、一个椭球体的名字。长半轴(semi-major axis)和反扁率(inverse flattening) 

5、本初子午线(prime meridian)名和其与格林威治子午线的偏移值 

6、投影方法类型（如横轴莫卡托） 

7、投影参数列表（如中央经线等） 

8、一个单位的名称和其到米和弧度单位的转换参数 

9、轴线的名称和顺序 

10、在预定义的权威坐标系中的编码（如EPSG） 

更多信息请参考OpenGIS WKT坐标系统定义，以及osr教程文档和OGRSpatialReference类的描述文档 

在GDAL中，返回坐标系统的函数是GDALDataset::GetProjectionRef()。它返回的坐标系统描述了地理参考坐标，暗含着仿射地理参考转换，这地理参考转换是由GDALDataset::GetGeoTransform()来返回。 

由GCPs地理参考坐标描述的坐标系统是由GDALDataset::GetGCPProjection()返回的。 

注意，返回的坐标系统字符串“”表示未知的地理参考坐标系统。 


2.1. 仿射地理转换
GDAL数据集有两种模式描述栅格位置（用点/线坐标系）以及地理参考坐标系之间的关系：首要的也是最普遍的是使用仿射转换，另一种则是GCPs（多控制点定位方式） 

仿射变换由6个参数构成，它们由GDALDataset::!GetGeoTransform()返回它们把点/线坐标系（我想这里的(点/线)意思是栅格的(行/列)位置）用下面的关系转换到地理参考空间。 


Toggle line numbers Toggle line numbers 
   1 Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
   2 Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

假设指北针向上的影像，GT2和GT4参数是0，而GT1是象元宽，GT5是象元高，（GT0，GT3）点位置是影像的左上角。 

注意，上面所说的点/线坐标系是从左上角(0,0)点到右下角，也就是坐标轴从左到右增长，从上到下增长的坐标系。点/线位置中心是(0.5,0.5) 


2.2. GCPs
数据集可以由一系列控制点来定义空间参考坐标系所有的GCPs共用一个地理参考坐标系(由GDALDataset::GetGCPProjection()返回，每个GCP（由GDAL_GCP定义）包含下面内容： 


Toggle line numbers Toggle line numbers 
   1  typedef struct
   2 {
   3     char        *pszId;
   4     char        *pszInfo;
   5     double      dfGCPPixel;
   6     double      dfGCPLine;
   7     double      dfGCPX;
   8     double      dfGCPY;
   9     double      dfGCPZ;
  10 } GDAL_GCP;

pszid是本GCP在数据集中一系列GCP点中惟一的标示字符串，它常常是数字，但不一定总是。有可能在GCP状态中包含机器可分析信息，虽然现在还不行。 

（dfGCPPixel，dfGCPLine）位置是栅格中的GCP位置,（dfGCPX，dfGCPY，dfGCPZ）位置是联合的地理参考位置（Z通常是0） 

GDAL数据模型没有实现由GCPs...产生坐标系的变化的机制，而是把它留给实际应用。但是1到5阶多项式是通常使用的方法。 

通常一个数据集会包含仿射地理变换。和GCPS中的一个或者两个都没有。两个都有很少见。而且无法用权威坐标系定义。 


3. 我的理解
上面文档说得很玄，还是让我们看看坐标系统究竟是什么样子。 


Toggle line numbers Toggle line numbers 
   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可以看这里 。 它是绝对的标准，完全解释了坐标系统的组成，结构和层次，它是由一些层次嵌套构成的。 


Toggle line numbers Toggle line numbers 
   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>}]