PostGIS之地理坐标

1. 概述

PostGIS 是PostgreSQL数据库一个空间数据库扩展，它添加了对地理对象的支持，允许在 SQL 中运行空间查询

PostGIS官网：About PostGIS | PostGIS

PostGIS官方教程：PostGIS 简介 — Introduction to PostGIS

PostGIS相关教程：文章目录汇总 - 知乎 (zhihu.com)

本文基于官方教程描述PostGIS中的地理坐标

数据准备可参考:

数据介绍可参考:

2. 地理坐标

2.1 平面坐标与球面坐标

地图中常用投影坐标系和地理坐标系，投影坐标系是将地球表面投影到平面上，地理坐标系是球面坐标，使用地理坐标系，不适用于运算面积、长度等

_images/cartesian_spherical.jpg

2.2 地理坐标

2.2.1 点到点

当使用的坐标系是地理坐标时，进行距离运算是没意义的，比如在WGS-84坐标系下计算洛杉矶（Los Angeles）和巴黎（Paris）的距离：

Los Angeles: POINT(-118.4079 33.9434)
Paris: POINT(2.3490 48.8533)

SELECT ST_Distance(ST_GeometryFromText('POINT(-118.4079 33.9434)', 4326), ST_GeomFromText('POINT(2.3490 48.8533)', 4326));

结果为：121.67388386839635，这似乎并没有什么意义

当使用的坐标系是地理坐标系时，需要将运算的背景放在球面坐标系，PostGIS中将这类数据定义为geography类型，而通常使用的平面坐标系下的属于geometry类型

使用geography类型来计算距离：

SELECT ST_Distance(ST_GeographyFromText('SRID=4326;POINT(-118.4079 33.9434)'), ST_GeogFromText('SRID=4326;POINT(2.3490 48.8533)'));

计算结果为：9125752.37394223，单位：米

2.2.2 点到线

在地理坐标系下进行距离量算，往往会出现错误的结果，比如下图中，"从洛杉矶（LAX）到巴黎（CDG）的航班路线距离冰岛有多远?"

使用平面坐标系，得到的是紫色的直线，而实际上，应该是红色的大圆路线

_images/lax_cdg.jpg

因此，这种情况下，应尽可能的使用geography类型来计算距离：

SELECT ST_Distance(
    ST_GeographyFromText('SRID=4326;LINESTRING(-118.4079 33.9434, 2.3490 48.8533)'),  -- LAX-CDG
    ST_GeogFromText('SRID=4326;POINT(-22.6056 63.9850)')  -- Iceland
);

2.2.3 国际日期变更线

国际日期变更线就是180°经线，根据定义，西经与东经相接于0°和180°经线，使用平面坐标系就会存在分界线问题，比如：WGS-84坐标系下，点（-179，0）与点（179，0）的距离是多少？

SELECT ST_Distance(ST_GeomFromText('SRID=4326;POINT(-179 0)'), ST_GeomFromText('SRID=4326;POINT(179 0)'));

结果为358°，显然这是错的，它计算的是（-179，0）到（179，0）的平面距离，没有考虑地理坐标的实际意义

因此，这种情况下，使用geography类型来计算距离：

SELECT ST_Distance(
    ST_GeogFromText('SRID=4326;POINT(-179 0)'), 
    ST_GeogFromText('SRID=4326;POINT(179 0)')
);

结果为：222638.98158655，单位：米，显然geography类型考虑了国际日期变更线的问题

2.3 Geometry转Geography

Geometry转Geography首先要将Geometry转换到地理坐标系WGS-84（SRID/EPSG：4326），再使用转换函数：

Geography(geometry) 将基于EPSG:4326的geometry数据类型转换为geography数据类型
（geometry）::geography 将基于EPSG:4326的geometry数据类型转换为geography数据类型

将Geometry数据表转换为Geography表：

CREATE TABLE nyc_subway_stations_geog AS
SELECT
  ST_Transform(geom,4326)::geography AS geog,
  name,
  routes
FROM nyc_subway_stations;

Geography表构建空间索引与Geometry相同，不同之处在于PostGIS对于Geography表会处理覆盖极点或国际日期变更线的要素的查询，而geometry空间索引则不会

CREATE INDEX nyc_subway_stations_geog_gix
ON nyc_subway_stations_geog USING GIST (geog);

对于Geography类型，支持的空间函数较少：

ST_AsText(geography) returns text
ST_GeographyFromText(text) returns geography
ST_AsBinary(geography) returns bytea
ST_GeogFromWKB(bytea) returns geography
ST_AsSVG(geography) returns text
ST_AsGML(geography) returns text
ST_AsKML(geography) returns text
ST_AsGeoJson(geography) returns text
ST_Distance(geography, geography) returns double
ST_DWithin(geography, geography, float8) returns boolean
ST_Area(geography) returns double
ST_Length(geography) returns double
ST_Covers(geography, geography) returns boolean
ST_CoveredBy(geography, geography) returns boolean
ST_Intersects(geography, geography) returns boolean
ST_Buffer(geography, float8) returns geography
ST_Intersection(geography, geography) returns geography

2.4 Geography转Geometry

与Geometry转Geography类似，支持转换的函数主要有：

Geometry(geography) 将基于EPSG:4326的geography数据类型转换为geometry数据类型
（geography）::geometry 将基于EPSG:4326的geography数据类型转换为geometry数据类型

SELECT Geometry(ST_GeogFromText('SRID=4326;POINT(-179 0)'));

2.5 创建Geography表

geography表可以创建时直接指定几何类型：

CREATE TABLE airports (
  code VARCHAR(3),
  geog GEOGRAPHY(Point)
);
 
INSERT INTO airports VALUES ('LAX', 'POINT(-118.4079 33.9434)');
INSERT INTO airports VALUES ('CDG', 'POINT(2.5559 49.0083)');
INSERT INTO airports VALUES ('KEF', 'POINT(-22.6056 63.9850)');

geography表的信息（元数据）储存在geography_columns表（视图）中