GIS开发扫盲贴--地理坐标系

GIS开发，可能最基础的，就是要搞懂坐标系了。

所谓GIS开发，无非就是处理一堆空间数据，即存储空间信息的数据。空间信息也者，经度纬度，高程之类。但是，地点就是那个地点，然后量度的方法、标准却多种多样，形成不同的所谓坐标系。

一、地理坐标系的分类

地理坐标系的分类与阐述，在百度百科里说得很复杂，阅读之后，只会让人脑袋嗡一声，傻半天。按照我目前粗浅的认识，我简单地把它归纳为两种：大地坐标系和投影坐标系。简而言之：

大地坐标系：
球面坐标系，以经纬度为单位。比如我们现在常说的某个地点，经纬度多少，就是这种球面坐标。GPS用的就是大地坐标系。

投影坐标系：

平面坐标系，常以米为单位。球面是曲面，但我们的地图是画在平面上的，这里要有一个换算过程。采取的就是投影法。比如墨卡托投影，用圆柱形的纸板将地球蒙起来，然后有一根蜡烛在地球球心，发光发热，将地球投射在这块纸板上，于是就得到了地球的平面图。

大地坐标系和投影坐标系是大类，下面又可以分为多种派系。比如投影坐标系，有墨卡托投影，高斯投影（高斯-克吕格投影），UTM投影。

二、WKID

这么多的坐标系，如何识别和区分呢？给它一个ID。在arcgis中，每种坐标系都有一个唯一标识，WKID。WKID，居然是Well Known ID，意为众所周知的ID。为什么说众所周知呢？因为这个ID不是arcgis独创的，它有国际标准，只不过，各个组织的叫法不同，但值是一样的。

1、OGC SRID
空间参考标识符 (SRID,Spatial Reference ID) 。

OGC(Open Geospatial Consortium)，开放地理空间信息联盟。自称是一个非盈利的国际标准组织，它制定了数据和服务的一系列标准，GIS厂商按照这个标准进行开发可保证空间数据的互操作。

OGC制订了许多标准：

Catalogue Service CS 用以发现、浏览服务器上数据、服务的元数据
CityGML 用以交换城市3D模型
KML KML 提供 XML 编码的地理数据集（从 Google 引入）
Coordinate Transformation Service CT 用以提供坐标系统及其转化的服务
Geography Markup Language GML 提供XML编码的地理数据集
Grid Coverage Service 栅格服务
Location Services LS 位置服务
Simple Features SFS 简单要素对象的通用描述
Web Coverage Processing Service WCPS 栅格处理Web服务
Web Coverage Service WCS 栅格Web服务 （Coverage栅格）
Web Feature Service WFS 要素Web服务（矢量）
Web Map Service WMS 地图Web服务
Web Map Tile Service WMTS 切片地图Web服务
Web Processing Service WPS 地理处理Web服务
Web Service Common OWS 描述了OGC Web服务的通用规范

其中有WKT(Well-known text)是一种文本标记语言，用于表示矢量几何对象、空间参照系统及空间参照系统之间的转换。它的二进制表示方式叫WKB(well-known binary)。

除了可以表示几何对象，还可以表示空间参照系统。一个表示空间参照系统的WKT字串描述了空间物体的测地基准、大地水准面、坐标系统及地图投影。 WKT在许多GIS程序中被广泛采用。ESRI亦在其shape文件格式(*.prj)中使用WKT。以下是空间参照系统的WKT表示样例：

COMPD_CS["OSGB36 / British National Grid + ODN",
    PROJCS["OSGB 1936 / British National Grid",
        GEOGCS["OSGB 1936",
            DATUM["OSGB_1936",
                spheroid["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],
                TOWGS84[375,-111,431,0,0,0,0],
                AUTHORITY["EPSG","6277"]],
            PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
            UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG","9108"]],
            AXIS["Lat",NORTH],
            AXIS["Long",EAST],
            AUTHORITY["EPSG","4277"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",49],
        PARAMETER["central_meridian",-2],
        PARAMETER["scale_factor",0.999601272],
        PARAMETER["false_easting",400000],
        PARAMETER["false_northing",-100000],
        UNIT["metre",1,AUTHORITY["EPSG","9001"]],
        AXIS["E",EAST],
        AXIS["N",NORTH],
        AUTHORITY["EPSG","27700"]],
    VERT_CS["Newlyn",
        VERT_DATUM["Ordnance Datum Newlyn",2005,AUTHORITY["EPSG","5101"]],
        UNIT["metre",1,AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","5701"]],
    AUTHORITY["EPSG","7405"]]

OGC用于标识不同空间参照系，用的是SRID，空间参考标识符 (SRID,Spatial Reference ID) ？

2、EPSG CODE
EPSG的英文全称是European Petroleum Survey Group，中文名称为欧洲石油调查组织。这个组织成立于1986年，2005年并入IOGP(International Association of Oil & Gas Producers)，中文名称为国际油气生产者协会。EPSG对世界的每一个地方都制定了地图，可能是为了找石油方便。

EPSG最初的目标是改善和共享世界各地的位置数据，负责发布并维护坐标参照系统的数据集参数，以及坐标转换描述，该数据集被广泛接受并使用。其中最著名的应该是EPSG:4326，美国主导的GPS系统就是在用它，它还有一个名气更大的别名叫作WGS84，WGS(World Geodetic System)是世界大地测量系统的意思，由于是1984年定义的，所以叫WGS84，之前的版本还有WGS72、WGS66、WGS60。

EPSG 用 EPSG CODE 来标识不同的坐标系。

3、ArcGIS WKID
Arcgis中通过WKID作为坐标参考系统的标识。

查找WKID实用网址：​​ArcGIS中的WKID​​

4、总结
SRID、EPSG CODE、WKID是不同组织和实体用于标识不同坐标系统的ID，但幸运的是，它们的值是一致的。

组织一多，可能就少不了磕碰、协调。

比如Web Mercator 坐标系使用的投影方法不是严格意义的墨卡托投影，而是一个被 EPSG称为伪墨卡托的投影方法，这个伪墨卡托投影方法的大名是 Popular Visualization Pseudo Mercator，PVPM。 看起来就觉得这个投影方法不是很严谨的样子，大众化的？受欢迎的？可视化伪墨卡托投影……因为这个坐标系统是 Google Map 最先使用的，或者更确切地说，是Google 最先发明的。在投影过程中，将表示地球的参考椭球体近似的作为正球体处理（正球体半径 R = 椭球体半长轴 a）。这也是为什么在 ArcGIS 中我们经常看到这个坐标系叫 WGS 1984 Web Mercator (Auxiliary Sphere)。Auxiliary Sphere 就是在告知你，这个坐标在投影过程中，将椭球体近似为正球体做投影变换，虽然基准面是WGS 1984 椭球面。后来，Web Mercator 在 Web 地图领域被广泛使用，这个坐标系就名声大噪。尽管这个坐标系由于精度问题一度不被GIS专业人士接受，但最终 EPSG 还是给了一个编号， WKID:3857。

而在此之前，arcgis已经给它分配了一个ID：102100（刚开始时是102113），那么现在就变成了在arcgis中两个ID都可以使用，都是指同一个坐标系。arcgis自己是这么标识的：

"spatialReference":{"wkid":102100,"latestWkid":3857}

三、坐标系转换

不同坐标系的空间数据，是不能混在一起工作的。比如，用arcgis for js叠加图层，两张图层的坐标系如果不一致的话，很可能是叠不上去的。除非arcgis默认会自动处理，比如wkid=4326与wkid=102100(3875)。

所以，使用之前，要将不同坐标系的数据进行转换，转成同一种。

怎么转换呢？以java为例：

import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.opengis.geometry.MismatchedDimensionException;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CRSAuthorityFactory;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;

public static Geometry trans2WGS1984(int wkid, Geometry geom)
        throws FactoryException, MismatchedDimensionException, TransformException {
        
    CoordinateReferenceSystem targetCRS = CRS.decode(String.format("EPSG:%d", 3857));
    /*
    强调经度在前。有些投影如高斯投影，XY与一般坐标系是相反的，
    X轴是我们通常意义上的Y轴，Y轴是X轴，刚好换过来。
    应对这种情况，就要使用这个函数。否则死活转换不成功，报什么：
    错误：The transform result may be 482.279 meters away from the expected position. 
    Are you sure that the input coordinates are inside this map projection area of validity? 
    The point is located 32°36.2'E away from the central meridian and 0°19.8'N away from the latitude of origin. 
    The projection is "Transverse_Mercator".
    */
    CRSAuthorityFactory factory = CRS.getAuthorityFactory(true);
    
    CoordinateReferenceSystem sourceCRS = factory.createCoordinateReferenceSystem(String.format("EPSG:%d",wkid));
    MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS,true);//最后一个参数为true有必要，表明是宽容模式，可以忽略高程
    
    return JTS.transform(geom,transform);
}

<dependency>
    <groupId>org.geotools</groupId>
    <artifactId>gt-shapefile</artifactId>
    <version>23.0</version>
</dependency>
<dependency>
    <groupId>org.geotools</groupId>
    <artifactId>gt-geojson</artifactId>
    <version>23.0</version>
</dependency>
<dependency>
    <groupId>org.geotools</groupId>
    <artifactId>gt-epsg-hsql</artifactId>
    <version>23.0</version>
</dependency>

也可以这样子来构造这个坐标系对象

= "PROJCS[\"CGCS2000 / Gauss-Kruger CM 123E\",GEOGCS[\"China Geodetic Coordinate System 2000\",DATUM[\"China_2000\",SPHEROID[\"CGCS2000\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"1024\"]],AUTHORITY[\"EPSG\",\"1043\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4490\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",123],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"4507\"],AXIS[\"X\",NORTH],AXIS[\"Y\",EAST]]";
  CoordinateReferenceSystem crsCGCS2000 = CRS.parseWKT(wktCgcs2000);

参考文章：
​​​arcgis for jsapi开发：坐标系、经纬度与平面坐标的互换​​​​Ogc标准​​​​WKT简介​​​​EPSG是什么？​​​​WGS_1984_Web_Mercator_Auxiliary_Sphere​​