GIS原理及应用、地理坐标系与投影坐标系

---

一、GIS定义

1.1 地理信息系统

1.1.1 简称GIS（Geographic Information System）。GIS四层模型如下。

• 表现层：直观展示地理数据，如数值化地图。
• 应用层：集合了对地理要素进行编辑、查询、分析的软件应用平台。
• 服务层：提供接口API，面向开发用户和最终用户提供数据处理、数据分析和发布服务。
• 数据层：体现了GIS的数据结构以及存储方式。

1.1.2 数据层
GIS数据库层，具备存储空间数据以及索引空间数据的能力。对应有文件型的ShapeFile、GeoJSON，ArcGIS的GeoDatabase，QGIS的GeoPackage，以及空间数据格式通过SDE-空间数据引擎架设到PostGIS/MySQL/Oracle上的相关数据等，通过空间数据引擎可以索引相关的数据。

1.1.3 服务层
服务层，提供接口API让二次开发人员、用户调用。如ArcGIS Server、QGIS Server、GeoServer等。

1.1.4 应用层
一些列GIS特有的功能和应用：空间数据加载、显示、编辑、查询、分析、处理和转换等。如ArcGIS Pro(ArcMap/ArcCatalog/ArcToolbox)、QGIS Desktop等。

1.1.5 表示层
如数字地图、行政区划图、专题图等。

1.2 建模

1.2.1 以排水行业为例进行说明，如下建模步骤。

• 抽象：将检查井抽象为点，管线抽象为线，泵站从整个管网研究尺度上抽象为点，再加上出水口等其他设施，将排水管网相关设施完成空间概化。
• 分层：虽然检查井和泵站都概化为点，但是由于各自管网模型中有不同角色和功用，且具有不同属性特征，故将其分为两类不同的要素，显示为不同的图层。
• 建模：构建排水管道专业模型-SWMM。GIS应用如下：
  • 流域分析（坡度分析），划分汇水区；
  • 网络分析，构建节点与管道的连接关系；
  • 叠加分析、提取高程信息。
    

1.3 相关教程

https://www.osgeo.cn/gis-tutorial/index.html
https://www.osgeo.cn/tutorial/
https://huzu.likegis.com/main/general/gis

1.4 GIS前沿方向

---

二、GIS数据格式

2.1 矢量

2.1.1 矢量数据就是代表地图图形的各离散点平面坐标 (x,y)的有序集合，主要用于表示地图图形元素几何数据之间及其属性数据之间的相互关系，并通过坐标的方式来记录点、线、面地理实体。

• 有序坐标集合（方向性）；
• 空间与属性关联；
• 坐标记录，点（0维）、线（1维）、面（2维）；

面的方向性：采用梯形法计算平面上任意三角形的面积，顺时针方向，如果采用终点的X-起点的X，那么面积为正；反之为负，如下：

2$S_{ABC}=(y_A+y_B)\ast (x_B-x_A)+(y_B+y_C)\ast (x_C-x_B)+(y_C+y_A)\ast (x_A-x_C)$
$=y_A{x}_B+y_B{x}_C+y_C{x}_A-y_B{x}_A-y_C{x}_B-y_A{x}_C$

2.1.1 可以采用GeoJSON、shp等格式存储；

2.2 栅格

2.2.1 栅格数据结构上就是像元的阵列，即像元按矩阵形式的集合，栅格中每个像元是栅格数据中最基本的信息存储单元，其坐标位置可以用行号和列号确定。

2.3 矢量与栅格的区别

• 连续性与离散性：栅格数据表示的地表是不连续的，是量化和近似离散的数据，这意味着地表一定面积内（像元地面分辨率范围内）地理数据的近似性。矢量数据其坐标空间为连续空间，不必像栅格数据结构那样进行量化处理，因此矢量数据更能精确地确定实体的空间位置。
• 属性值与像元值：矢量数据一般是现实世界实体的图形化抽象，因此其关联的属性数据为对应的抽象实体的属性信息；栅格数据的属性值一般对应波段的像元值；
• 图形学中的矢量数据，如svg，而GIS中的矢量数据可近似理解为图形学矢量数据+属性信息+投影信息；图形学中的栅格数据，如BMP，而GIS中的栅格数据可近似理解为图形学栅格数据+投影信息。

---

三、GIS数据组织

3.1 抽象

3.1.1 空间层面，根据研究尺度和应用需求按实际事务抽象为点、线、面要素。

3.1.2 属性层面，数据结构设计，抽象要素（Feature）的属性；

3.2 分层

3.2.1 图层（矢量图层Vector Layer）（FeatureLayer）：将具有同类特征的要素（Feature）归为一个要素类（FeatureClass），通过一定的配置（标注、渲染符号等）渲染到地图。

3.2.2 图层的组成与分解

• 矢量图层（FeatureLayer）
  • 数据源要素类（FeatureClass）
    • 元数据字段结构Field；
    • 数据表要素Feature
  • 图层装饰
    • 图层标注Label
    • 数据时态Timing
    • 符号渲染Symbology
    • 提示框Tooltip
    • 数据过滤Filter
    • 数据关联
    • …

3.2.3 矢量图层数据源格式主要有：

• ESRI ShapeFile（.shp）标准格式；
• GeoJON常见于WebGIS数据交换；
• ESRI GeoDatabase
• QGIS GeoPackage

3.3 栅格与切片

3.3.1 遥感影像
空间分辨率（Spatial Resolution）又称地面分辨率；

• 后者是针对地面而言，指可以识别的最小地面距离或最小目标物大小；
• 前者是针对遥感器或者图像而言的，指图像上能够详细区分的最小单元的尺寸或大小；

卫星遥感图像的计量单位为景，即一景指卫星拍摄一次所成像的画面，一景的范围一般在几千至上万平方公里。

3.3.2 栅格数据除了常见的遥感影像，还有一类叫数字地形模型DTM，其中一个分支就是DEM数字高程模型。

• 数字高程模型DEM：
  • 线表示：等高线；
  • 点表示
    • 侠义的DEM，如规则网格Grid
    • 不规则三角网TIN；

3.3.3 切片

地图切片（地图瓦片），瓦片数据是将矢量数据或影像数据进行预处理，采用高效的缓存机制（如金字塔）形成缓存图片集，采用“级、行、列”方式进行组织，可在网页中快速加载。

3.3.4 地图切片

• 切片原点：左上角、中点；

  • 大部分切片算法是采用左上角为切片原点，如天地图、高德地图、谷歌地图等。切片原点定义了第一个256*256的切片的位置。

• 切片分辨率（比例尺）

  • 切片分辨率和比例尺是一一对应的。比例尺是指地图上一个单位代表实际多少距离，而分辨率代表屏幕上一个像素代表实际多少距离，因此，当确定一个后，对应的另一个也就确定了。

• 单个切片像素大小：一般采用256*256px；

---

四、坐标系

4.1 坐标系简介

4.1.1 EPSG，即European Petroleum Survey Group，欧洲石油调查组织，它为每个坐标系都分配了一个介于1024~32767的数字代码。

4.1.2 Web Mercator，EPSG:3857,即Web墨卡托投影，投影后整幅世界地图成正方形。不同的坐标投影出来的地图形状各异。

4.1.3 不同的坐标系有不同的特点，有些投影能保持角度不变，但距离和面积上会有偏差，比如纬度越高的地方偏差越大。

4.2 大地坐标系GCS

4.2.1 大地坐标系。也称地理坐标系（GCS）。

4.2.2 大地水准面：假定海水处于“完全”静止状态，把海水面延伸到大陆上去，形成包围整个地球的连续表面。
4.2.3 地球椭球面：凡与局部地区（一个或几个国家）的大地水准面符合的最好的旋转椭球，称为“参考椭球”。

4.2.4 地心坐标系与参心坐标系

• 地心坐标系：选用地球质心作为椭球中心，从拟合全球曲面的角度出发来建立空间坐标系，适用于全球，如WGS84、CGCS2000。
• 参心坐标系：从拟合局部区域的角度出发来建立空间坐标系，选用最贴合局部地区曲面的椭球中心，适用于地方，如西安80、北京54。
  

4.2.5 地球椭球面（Datum）

• 椭球体参数
  • 长半轴a
  • 短半轴b
  • 扁率alpha
    $alpha=\frac{(a-b)}{a}$
  • 第一偏心率e
    $e=\sqrt{\frac{a^2-b^2}{a^2}}$
  • 第二偏心率e‘
    $e^{\prime }=\sqrt{\frac{a^2-b^2}{b^2}}$
    
    WGS84椭球体的参数为：a=6378137，alpha=$\frac{1}{298.257223563}$,CGCS2000椭球体的参数为：a=6378137，alpha=$\frac{1}{298.257222101}$,根据参数a、alpha就可计算其他参数。

4.2.6 大地坐标系，有两种表示法，分别为采用大地坐标经纬高：（L，B，H）和空间直角坐标系：（X，Y，Z）来表示。

• P点向球体做法线，可知（X，Y，Z）与（L，B，H）之间的关系，由于是椭球体，长半轴与短半轴不同，因此法线的延长线与空间直角坐标系不会相交于原点（球心）。

• 在不考虑小数精度的情况下，此二者转换不会丢失精度。

• 在地球表面时，H为0；

• 零纬度线称为赤道；

• 零经度线称为本初子午线。对于绝大多数地理坐标系，本初子午线是指通过英国格林尼治的经线；

• 经纬网的原点（0，0）定义在赤道和本初子午线的交点处。

4.3 投影坐标系PCS

4.3.1 投影坐标系=大地坐标系+投影方法（投影函数或映射函数）；

• 平面坐标（X，Y）= F（L，B），F为函数映射，L为经度、B为维度；
  
  为了便于理解几何透视原理，想象地球是一个表面透明的球体，其上绘有经纬网，用一张巨大的纸（称为投影曲面）包裹地球，假设有一个位于地心处的光源穿过地球将经纬网投影到这张上，然后将其展开就可以得到一张地图。

4.3.2 投影坐标系是建立在一种大地坐标系之上的。

4.3.3 投影坐标系分类

• 按照变形性质分为:等角投影、等距、等积投影等；
• 按照地球与投影面的相对位置称为：横轴、正轴、斜轴投影；
• 按照投影面分为：圆柱、圆锥、平面投影等；
  

4.3.4 Web Mercator（EPSG：3857）墨卡托投影

• 大地坐标GCS为：WGS84（EPSG：4326）；
• 投影方法：采用正轴等角圆柱投影，即墨卡托投影。
  • 等角；
  • 正轴；
  • 圆柱投影面；

Web Mercator = WGS1984 + 正轴等角圆柱投影，投影结果为正方形。

投影误区：
投影并非直线照射投影到柱面上，即非简单的几何透视投影。采用反证法可以证明，已知当B约等于85.05度，投影展开为正方形，如简单的理解为透视投影，则正方形的长为赤道长度，宽为2L，有：
$长=2\ast PI\ast R，宽=2L=2tanB\ast R，则：B=72.343度$
B与85.05度不相等，因此不能简单的理解为几何透视投影。

4.3.5 CGCS2000 Gauss-Kruger CM 117E(EPSG:4509)

• GCS:CGCS2000（EPSG：4490）
• 投影方法：等角横切圆柱投影，即高斯-克吕格投影；
  • 等角；
  • 横轴；
  • 圆柱投影面；

CGCS2000 Gauss-Kruger CM 117E= CGCS2000 + 高斯-克吕格投影（6度带，以东经117为中央子午线），即117度加减3度；

4.3.6 投影比较

	Web Mercator	CGCS2000 CK CM 117E
大地坐标系	WGS1984	CGCS2000
投影方式	墨卡托（等角正轴切圆柱）	高斯-克吕格（等角横轴切圆柱）
使用场景（优势）	等角（航海）；适用于小比例尺，全世界地图；适合Web，利于地区图切片	等角（城市交通、道路交叉）；适用于大比例尺，城市级别；高斯克吕格投影在长度和面积上变形很小，中央经线无变形
劣势	纬度越高，变形严重；面积无参考价值；	不适合大范围，全国地图涉及拼接，而每带为瓜瓣状，难以无缝拼接

4.4 投影变换

4.4.1 通过分析源代码，理解投影变换。
https://github.com/proj4js/proj4js，transform.js源代码如下。

import {D2R, R2D, PJD_3PARAM, PJD_7PARAM, PJD_GRIDSHIFT} from './constants/values';
import datum_transform from './datum_transform';
import adjust_axis from './adjust_axis';
import proj from './Proj';
import toPoint from './common/toPoint';
import checkSanity from './checkSanity';function checkNotWGS(source, dest) {return ((source.datum.datum_type === PJD_3PARAM || source.datum.datum_type === PJD_7PARAM || source.datum.datum_type === PJD_GRIDSHIFT) && dest.datumCode !== 'WGS84') ||((dest.datum.datum_type === PJD_3PARAM || dest.datum.datum_type === PJD_7PARAM || dest.datum.datum_type === PJD_GRIDSHIFT) && source.datumCode !== 'WGS84');
}export default function transform(source, dest, point, enforceAxis) {var wgs84;if (Array.isArray(point)) {point = toPoint(point);} else {// Clone the point object so inputs don't get modifiedpoint = {x: point.x,y: point.y,z: point.z,m: point.m};}var hasZ = point.z !== undefined;checkSanity(point);// Workaround for datum shifts towgs84, if either source or destination projection is not wgs84if (source.datum && dest.datum && checkNotWGS(source, dest)) {wgs84 = new proj('WGS84');point = transform(source, wgs84, point, enforceAxis);source = wgs84;}// DGR, 2010/11/12if (enforceAxis && source.axis !== 'enu') {point = adjust_axis(source, false, point);}// Transform source points to long/lat, if they aren't already.if (source.projName === 'longlat') {point = {x: point.x * D2R,y: point.y * D2R,z: point.z || 0};} else {if (source.to_meter) {point = {x: point.x * source.to_meter,y: point.y * source.to_meter,z: point.z || 0};}point = source.inverse(point); // Convert Cartesian to longlatif (!point) {return;}}// Adjust for the prime meridian if necessaryif (source.from_greenwich) {point.x += source.from_greenwich;}// Convert datums if needed, and if possible.point = datum_transform(source.datum, dest.datum, point);if (!point) {return;}// Adjust for the prime meridian if necessaryif (dest.from_greenwich) {point = {x: point.x - dest.from_greenwich,y: point.y,z: point.z || 0};}if (dest.projName === 'longlat') {// convert radians to decimal degreespoint = {x: point.x * R2D,y: point.y * R2D,z: point.z || 0};} else { // else projectpoint = dest.forward(point);if (dest.to_meter) {point = {x: point.x / dest.to_meter,y: point.y / dest.to_meter,z: point.z || 0};}}// DGR, 2010/11/12if (enforceAxis && dest.axis !== 'enu') {return adjust_axis(dest, true, point);}if (point && !hasZ) {delete point.z;}return point;
}

代码流程图如下：

• 确保源坐标系或目标坐标系有一个是WGS84，若都不是，则先把目标坐标系转换为WGS84，在从WGS84转换到目标坐标系；
• 如果源坐标是投影坐标系，则将平面坐标转换为经纬度坐标；
• 如果源坐标系和目标坐标系椭球体（基准面）参数不一样，则需要进行椭球体参数之间的转换；这两个椭球体对应的空间直角坐标系之间，可以通过三参数或七参数转换方法来进行转换；

两个不同的大地坐标下不同的投影坐标系下的转换，简化后的变换流程大致如下：

4.4.2 空间直角坐标系-7参数或三参数转换

• 七参数转换：三个平移参数、三个旋转参数和一个缩放参数M，转换精度更高，适用于区域范围更广，控制点起码3个以上。
  
• 三参数转换：三个平移参数，适用范围小，最少需要一个控制点。
  

4.4.3 Web Mercator变迁史

• 别名：伪墨卡托（Pseudo Mercator）
• WKID：3857；
• 基准面：WGS1984；
• 投影方法：正轴等角圆柱投影；
• 投影范围：经度$-180^{。}$~$180^{。}$，南北纬$85.05^{。}$
• 争议：采用椭球体的坐标系，却用球体的参数来投影计算；导致问题：并非等角，高纬度问题更显著。
  

4.4.4 GCJ02简介
主要用于对在线发布的地图进行坐标的偏移，采用了GCJ02偏移算法，偏移的是经纬度，而不是平面坐标，GCJ02不应该称为坐标系，而应该称之为偏移算法，GCJ02并没有基于任何大地坐标系进行变换。

---

五、空间数据库与空间数据索引

5.1 空间索引

5.1.1 空间索引是指空间要素的地理位置，形状或空间对象之间的某种空间光系，按一定的顺序排列的一种数据结构，一般包括空间要素标识，外包络矩形以及指向空间要素的指针。

5.1.2 外包络矩形（Envelope）
外包络矩形是指空间要素的封装边界。

5.1.3 矩形范围索引
基本原理就是对空间要素的外包络矩形进行索引。在在进行空间查询时，根据空间要素外包络矩形与查询范围进行比较，判断是否有叠盖而得到查询结果。

5.1.4 单位网格索引
基本原理就是对空间要素的外包络矩形进行索引，在进行空间范围查询时，分为两级过滤（筛选）。初次过滤根据空间要素外包络矩形来过滤掉大部查询的空间要素，因为空间要素外包络矩形已被索引，所以初次过滤过程比较快，花费的代价较小。第二级过滤则用查询空间范围直接和初次过滤结果集中空间要素的二进制边界坐标比较，从而查询到准确的结果。

5.1.5 单位网格索引改进
采用二维格网索引+外包络矩形来进行空间索引。

5.1.6 R树索引

B树的空间版，根据空间位置来维护插入和删除的B树。

5.1.7 四叉树索引

5.1.8 总结
https://www.osgeo.cn/gis-tutorial/ch06-06/sec05-5.html

六、经验之谈

6.1 在GIS中，投影坐标系下模型顶点坐标[x,y]值的范围及数量级通常与具体的投影方法和分带有关。在中国范围内，通常采用高斯-克吕格投影，其中x坐标通常是一个较大的值，而y坐标则是一个相对较小的值。y坐标通常是7位的，而x坐标有两种情况，6位和8位。6位x坐标通常省略带号，是中央经线形式的投影坐标，而8位x坐标带有带号。

6.2 在高斯投影直角坐标中，为了避免出现负的横坐标，横坐标一般会加上500km，意味着x坐标通常是以十万（500000米）位单位的，而y坐标的最大值（在赤道上，6度带）约为330公里。

附录

https://www.bilibili.com/video/BV1Sa411V7aM?spm_id_from=333.788.videopod.sections&vd_source=d5f54bdbe1c20d376b8816eaa36c75b0