浅谈geometry和geography
因为论文中使用postgis作为空间数据库拓展，然后自己有点搞混这两个东西，所以在这里提一嘴。
当然有些内容可能有所偏差，并不能保证正确性，如有异议，可以提出来。
起因
起因是有一个需求，数据是经纬度，用geometry存储，找出近多少米范围内的所有点（或者线面），所以就去找一个函数，找到这篇教程PostGIS教程九：空间关系（顺带一提，这个系列教程对于初学者使用postgis很有用），里面介绍了一个ST_DWithin的函数，用了一下发现不对，就随便设个半径都能有大量的数据，查了一下官方文档ST_DWithin，发现是单位的问题，在此坐标系中，用geometry是以度为单位，而两点之间的距离用度为单位没有什么意义，比如极限一点，赤道的一度和极点的一度不是一个概念，所以需要转换成geography，这样子的单位就是米了，才能正常使用。但是我对于这块的理解还不是很到位，一开始不是很在意，在论文里提了一嘴，但是老师说我写的不对，最后争论半天，还是把这块内容给去除了，最近小论文实验又需要调整，所以这块内容再次省视。
坐标系
这块内容直接抄的PostGIS中geometry与geography的区别这篇博文，讲的非常好。
同时也参考这篇文章blog/20171018_02.md，这个带佬是postgresql方面的带佬，讲了很多篇postgresql。
看图讲话，就是GIS有两种坐标系（不确定是不是只有这两种），一种是平面坐标（投影坐标、笛卡尔坐标），就是下图左边这种，是二维平面使用的，以前数学的几何题目都是这种坐标系类型。另一种是球坐标（大地坐标、地理坐标、经纬度坐标），就是下图右边这种，经纬度采用的就是这种坐标系。
投影坐标就是球坐标投影展开后得来（就是圆柱把🌏包起来，然后地球是发光体，投影后再把圆柱展开）
目前用得最多的有SRID=4326球坐标（另一个大气的别名就是WGS84，就是美国GPS使用的坐标系，也是目前使用最为常见和广泛的），SRID为EPSG:3857的墨卡托投影坐标（就是上述的投影之后的二维坐标系）。
我国自主的大地坐标系是EPSG:4490，单位也是经纬度，另一个别名是CGCS2000，它是最适合中国地理位置的地理坐标。当然以前是北京54地理坐标和西安94地理坐标，不过已经逐渐变成历史了，现在就是2000国家大地坐标系。
这个网站能查各种坐标系Spatial Reference以及这个http://epsg.io/
也可以在postgis中查询与北京、中国有关的坐标系。
select * from spatial_ref_sys where srtext ~* 'beijing';
select * from spatial_ref_sys where srtext ~* 'china';
比如我想把EPSG:4326的地图映射到比较符合北京市的平面投影坐标系，就可以去找，北京在东经115-117度，所以我找到EPSG:4796这个比较合适，这是我国新的平面投影坐标，如下图所示。之前的北京54类似的投影坐标系 EPSG:2436，是有纬度的限制，但是4796没有纬度限制，我简单测试了一下这两个，发现没有区别，所以有新的就用新的投影坐标系吧。
几何和地理
geometry翻译过来应该就是几何，geography就是地理。一开始我以为就是分别对应几何坐标和球坐标，但是不完全是。geometry不仅可以使用几何坐标，同时也可以使用球坐标。
这里用抄来的表格做个总结。
| 对比项 | geometry | geography |
|---|---|---|
| 名称 | 几何对象 | 地理对象 |
| 坐标系 | 支持平面坐标系和球面坐标系 | 仅支持球面坐标系 |
| 单位 | 视使用的坐标系的情况，不同的坐标系的单位是不一样的，比如EPSG:3857使用的就是米为单位，但是EPSG:4326的单位是度 | 单位是度？（不确定，因为有其他球坐标系，极坐标？，可能单位不是度，但是一般都是经纬度，所以单位是度比较默认常见） |
| 对象类型 | 支持 POINT、MULTIPOINT、LINESTRING、LINEARRING、MULTILINESTRING、POLYGON、MULTIPOLYGON、POLYHEDRALSURFACE、TRIANGLE、TIN、GEOMETRYCOLLECTION等简单对象，还支持CIRCULARSTRING、COMPOUNDCURVE、CURVEPOLYGON、MULTICURVE、MULTISURFACE | 仅支持POINT、LINESTRING、POLYGON、MULTIPOINT、MULTILINESTRING、MULTIPOLYGON、GEOMETRYCOLLECTION |
| 函数限制 | 类型较丰富 | 支持类型较少，仅支持类型转换、长度面积距离计算、交并差运算函数。可以参考这篇PostGIS教程十三：地理 |
| 索引局限性 | 几何索引不能正确处理极地地区查询（毕竟是投影过来的，上面那块相当于没有投影过来，所以不能正确处理） | 地理索引可以正确处理覆盖极点或国际日期变更线的查询 |
| 元数据 | geometry_columns：提供了数据库中所有空间数据表的描述信息，分别是数据库名、模式名、空间数据表名，属性列名称，几何图形维度，空间参考标识符，几何图形类型 | 无定义的元数据，需自定义 |
从上面的表格看出来，geometry更加强大，不仅能够支持几何坐标，也能支持球坐标，而且空间函数更加丰富。
使用选择可以参考以下表格
| geometry | geography |
|---|---|
| 如果数据在地理范围上是紧凑的（包含在州、县或市内），推荐使用基于笛卡尔坐标的geometry类型（因为比较紧凑，投影的失真更少，所以误差会更小） | 如果需要测量在地理范围上是分散的数据集（覆盖世界大部分地区）距离，推荐使用geography类型。（比如欧洲和美洲，虽然geography计算比较费时，但是但是这么远的距离能够抵消这种计算的费时） |
| 当做数据存储时，推荐使用geometry | 由于地理坐标较为精确，因此在进行距离、面积等量算时，建议使用geography（注意：经纬度两者之间的距离计算一定要换到地理坐标！！） |
| 如果用户较为了解投影信息知识，推荐使用geometry | geography不需要了解专业的投影知识只需要知道经纬度就可以进行计算，因此使用门槛较低（就是经纬度，谁都会用） |
| 当场景需要运用大量复杂函数时，推荐使用geometry | geography的支持函数较少，且计算复杂，因此应用时需要占用较多计算资源。 |
从上述看出来，两者的差距还是有点的，尤其是在进行距离计算的时候，必须是要转换成geography。
互相转换
两者之间可以相互转换，因为geometry支持两个坐标系，所以先需要将坐标系转换到球坐标系，也就是EPSG:4326（经度/纬度），然后再把格式转换成geography，就如下图所示。当然如果你的geometry存储的就是就是4326的坐标，那么只需要第二步转换就行了。
那么逆向的转换也是一样的，反过来就行了。
更多的转换和使用姿势可以看这篇如何理解类型geometry和geography以及4326、3857坐标系。
正确使用
1、geometry 球坐标
对于geometry类型，可以使用st_distancespheroid直接用球坐标，第三个参数通过指定SPHEROID来精确计算（这块不是很懂，大概就是球体的质心和轴长每个地方有点差异，可以输入当地比较准的）。
SELECT st_distancespheroid(ST_GeomFromText('POINT(116.34 39.92)', 4326),
ST_GeomFromText('POINT(116.34 39.96)', 4326), 'SPHEROID["WGS84",6378137,298.257223563]');
2、geometry 平面坐标
也就是先把球坐标投影到比较精准的投影坐标（小面积投影坐标系），但是必须要覆盖到两个点。
我这里给的两个点都在北京市，所以用ESPG4796就能够覆盖这两个点，用这个比较符合北京的平面坐标系误差就会小，然后使用普通的平面求距离的st_distance就可以了。
提示：如果不是跨度非常大的坐标，比如跨洲啊这种，那么建议转换成符合坐标情况的平面坐标系，单位为米的那种，因为geometry的相关操作很快，比起geography的三角计算快的太多了。
select st_distance(ST_Transform(ST_GeomFromText('POINT(116.34 39.92)', 4326), 4796 ),
ST_Transform(ST_GeomFromText('POINT(116.34 39.96)', 4326), 4796 ));
3、geography 球坐标系
这里显然是球坐标系，应该也用不了其他的坐标系，但是球坐标系好像也有好多种，除了4326，但是问题不大，对于ST_Distance的第三个参数，可以使用true，那么就可以适用于椭球体（WGS84）
float ST_Distance(geography gg1, geography gg2, boolean use_spheroid);
注意这里不是ST_GeomFromText，而是ST_GeogFromText
SELECT st_distance(ST_GeogFromText('SRID=4326;POINT(116.34 39.92)'),
ST_GeogFromText('SRID=4326;POINT(116.34 39.96)'), true);
解决问题
方法一
回到最开始我的需求，找出某个圆内的所有点，应该怎么弄呢
首先是对这个数据库启用postgis扩展
create extension postgis;
然后创建表，这里采用的是geometry格式，因为更加通用（如果都是经纬度操作的话，用geography格式其实更合适）
create table t_pos(
id int primary key,
geom geometry(geometry,4326)
);
创建索引，注意这里创建了两条索引，一个是正常的geometry的索引，另一个是geography的索引，如果不创建geography的索引，那么使用ST_DWithin需要转换成geography的时候，他会遍历所有表中的数据转换成geography，这样子速度会很慢，所以需要建立geography索引，虽然在添加的时候会慢一点。
create index gidx_t_pos_geom on t_pos using gist(geom);
create index gidx_t_pos_geog on t_pos using gist(geom::geography);
然后使用ST_DWithin来查找附近5000米以内的所有数据，注意这里函数的第一个和第二个参数都需要是geography类型的才行，可以使用::geography进行转换，函数的第三个参数use_spheroid为true表示使用椭球体（spheroid为椭球的意思），如果为false的话，就是使用球体计算（sphere为球体），精度下降，但速度更快，如果一二参数都是geography ，那么这个值默认为true，不写也没事。
select * from t_pos where ST_DWithin(geom::geography,ST_SetSRID(ST_Point(102,24),4326)::geography,5000,true)
方法二（推荐）
同样开启postigs拓展，然后建表，这里建表就不一样了，这里指定里面的坐标系是平面坐标系ESPG:4796，然后创建一条索引就行了。
create table t_pos(
id int primary key,
geom geometry(geometry,4796)
);
create index gidx_t_pos_geom on t_pos using gist(geom);
当然插入数据的时候也需要转换相应的坐标系，这里从4326球坐标转换到4796平面坐标系
insert into t_pos values (1, ST_Transform(ST_GeomFromText('SRID=4326;POINT(116.34 39.92)'), 4796));
insert into t_pos values (2, ST_Transform(ST_GeomFromText('SRID=4326;POINT(116.34 39.96)'), 4796));
然后测一下距离先，4441.5899530316865是正确的
SELECT st_distance((select geom from t_pos where id=1),(select geom from t_pos where id=2));
然后我们就可以愉快的用ST_DWithin，而不需要转换成为geography。
116.34经度下，北纬39.95离39.96为1110米，离39.92为3331米，所以可以尝试1112和3330和3332这几个值，看是否符合实际预期。
select * from t_pos where ST_DWithin(geom,
ST_Transform(ST_GeomFromText('SRID=4326;POINT(116.34 39.95)'), 4796),1112,true)
疑问
既然geometry支持4326球坐标，为什么还要geography来计算距离呢，难道不能直接使用geometry吗？这块我没有去仔细查，我的理解可能应该是虽然geometry支持这个4326，而且是以度为单位，但他的处理方法还是以二维平面的方式去处理的，并且单位还是度，所以在ST_DWithin这种函数中，我们如果参数为geometry，他会以geometry坐标系给的单位进行计算，给啥就是啥，如果是度就用度（因为他不知道你用的是平面坐标还是球坐标吧？）。而是用geography就一定是球坐标，所以他会使用geography的方法进行计算，是用的函数计算都是通过大量的三角运算，对于距离这种，就会转换为用米来作为单位进行计算了。
但是对于KNN这种来说，其实应该也用geography来计算的，但是我目前说实话，没有看到有人这样子使用，都是用的geometry来计算knn，我以前测试过，使用geography来计算knn速度会非常的慢，但是用geometry就非常得快，当然我没有看具体的内容，所以到底geometry的球坐标使用knn是否是符合实际的，他们内部是否有做过优化，或者说用度作为单位来knn会有很大的影响吗，这点目前还不清楚，明天尝试一下。
尝试了，geometry的knn确实不行，应该是按照度为单位来计算的，所以并不是真正的距离从近到远，但是粗略的使用也还是可以吧。
所以最佳实践是什么呢？就是将ESPG:4326球坐标系映射到符合数据情况的平面坐标系，又快又好又稳。
