基于PostGIS的火星坐标，百度坐标批量互转
一 背景
GIS项目中底图是必不可少的，绝大部分GIS项目使用的底图是基于高德，谷歌，百度，天地图等互联网（在线/离线）底图。由于我国特殊国情，公众版地理信息服务（包括电子底图）都要进行各种坐标偏移旋转等数据加密处理，并获取国土资源部数据审查并颁发审图号才可公开发布。企业的GIS数据都是通过传感器或者实地测量获得的非加密的WGS84坐标（即常用的gps那种坐标），当企业将自己的业务数据叠加到互联网底图时，不可避免出现图层叠加偏移问题，如下图：
图层叠加偏移.png
底图偏移情况总结如下：
- 百度底图，最坑爹的，bd-09坐标，二次加偏。
- 高德，谷歌底图，GCJ-02坐标，俗称的“火星坐标系”。
- 天地图，osm底图，通常认为是不偏移底图。
二 解决方案与存在问题
2.1 解决方案
由于底图是不可改变的，唯一能解决叠加偏移问题的方法是：将自己的业务数据根据一定的偏移公式转换，从而使数据的坐标与底图达到相对对齐。当前主流编程语言，都有类似处理坐标偏移的公式换算库，满足不同语言开发的GIS项目需求。比如如下几个地址：
https://github.com/FreeGIS/CoordinateTransform
https://github.com/wandergis/coordtransform
2.2 存在问题
- 不支持批量计算。
当前各种坐标转换库，都是实现单个点坐标转到底图对应的坐标。但实际上，GIS存在 图形复杂，数据量大的通用特征，每个图层都有大量图形记录，每个图形记录，其图形常用分类就有点线面多点多线多面等六种，每个图形都可能由非常多的坐标点组成。业务需求最好是支持批量计算。 - 动态计算，显示效率低。
当前各种坐标转换库，都是基于前后台的语言实现的，每次都要获取数据源数据，再通过实时计算展示，当数据复杂，数据量大时候，效率很成问题。解决办法通常是，后台先读取数据源，然后计算好，最后将计算结果重新存入库里，项目使用时，直接读取计算好的结果。 捯饬的费劲不？ - 数据转换服务维护复杂。
当数据库数据发生变更时，转换服务要能增量得到变更的数据，然后计算，然后再更新数据库已存储的转换结果。整个架构要是能自动维护好，还是很复杂的。
三 PostGIS方案
基于现有方案普遍存在的问题，笔者基于最常用的开源GIS空间数据库PostGIS开发了一个function，基于一行sql搞定一切的是思想，满足使用PostGIS的开发者，简单的处理这些坐标转换问题。PostGIS实现的仓库如下：
https://github.com/FreeGIS/Postgis_Coordinate_Transform
实现功能列表如下：
- 支持WGS84与bd-09,gcj-02坐标系，百度经纬度与百度墨卡托之间互转。
- 支持点线面多点多线多面的复杂图形批量转换。
- 支持对整个表批处理转换。
约束：
- 要求转换的表是基于PostGIS创建的空间关系表
示例支持的表：
create table point_test(
gid serial primary key,
name text,
geom geometry(Point,4326)
);
不支持的表：
create table point_test(
gid serial primary key,
name text,
lon numeric, --经度
lat numeric --纬度
);
不支持的表是普通关系表，非空间图形表。
- 要求转换的图形必须是二维图形
当前暂不支持三维或者多维，如带Z值的高程，带M值的测量值等，由于过于复杂，目前笔者暂未有时间去过多实现。 - 转换表图形坐标系必须是epsg:4326
除了将百度墨卡托坐标转百度经纬度外，其他转换方式，必须保证转换表的坐标系是4326，其他坐标系，需要用户使用ST_Transform函数，将其数据先转到4326坐标系下，再使用该工具。百度墨卡托坐标转百度经纬度转换，数据源必须是3857的。 - 图形数据是点线面多点多线多面
仅仅支持Point,LineString,Polygon,MultiPoint,MultiLineString,MultiPolygon六种明确类型。其他的PostGIS类型由于不常用，且不严格规范，通常不用于标准的空间数据库类型，暂时不考虑实现。
3.1 安装应用
前提：PostGIS用户，图形表是基于PostGIS的空间关系表。
示例：在test库安装转换方法
[postgres@sss~]$ psql -d test
psql (11.1)
Type "help" for help.
test=# \i FreeGIS_Coordinate_Transform.sql
CREATE FUNCTION
CREATE FUNCTION
CREATE FUNCTION
CREATE FUNCTION
CREATE FUNCTION
CREATE FUNCTION
CREATE FUNCTION
CREATE FUNCTION
CREATE TYPE
CREATE FUNCTION
方法简介：
FreeGIS_Coordinate_Transform(
in schema_name text,
in table_name text,
in transform_type FreeGIS_coordinate_transform_type
)
schema_name:表的schema名称。
table_name：表的名称。
transform_type：转换类型，枚举型，类型如下：
BD2GCJ：百度经纬度 转 火星经纬度。
GCJ2BD：火星经纬度 转 百度经纬度。
WGS2GCJ：WGS84经纬度 转 火星经纬度
GCJ2WGS：火星经纬度 转 WGS84经纬度。
BD2WGS：百度经纬度 转 WGS84经纬度。
WGS2BD：WGS84经纬度 转 百度经纬度。
BDWGS2BDMKT：百度经纬度 转 百度墨卡托。
BDMKT2BDWGS：百度墨卡托 转 百度经纬度。
WGS2BDMKT： WGS84经纬度 转 百度墨卡托。
BDMKT2WGS：百度墨卡托 转 WGS84经纬度。
转换批处理执行：
--将public.test表从wgs84坐标转火星坐标。
select FreeGIS_Coordinate_Transform('public','test','WGS2GCJ');
执行后，test表新增了一个transform_geom字段，就是转换后的结果。
转换结果
3.2 底图叠加示例
3.2.1 叠加osm，天地图底图
企业的坐标，可不做处理，直接叠加。
点叠加osm
3.2.2 叠加谷歌，高德底图
将企业的数据，使用wgs84转火星坐标，然后叠加：
select FreeGIS_Coordinate_Transform('public','test_pt','WGS2GCJ');
执行展示：
--transform_geom字段存储了转换结果，示例使用如下：
create or replace view v_test_pt as select gid,name,transform_geom as geom from test_pt;
将v_test_pt 发布成地理服务后，叠加地图如下：
转换后叠加osm，发生偏移
转换后叠加谷歌或高德，位置匹配
依次制作了测试线和面，测试情况如下
线 osm
线 谷歌
面 osm
面 谷歌
3.2.3 叠加百度底图
将企业的数据，使用wgs84转百度墨卡托坐标，然后叠加：
select FreeGIS_Coordinate_Transform('public','test_pg','WGS2BDMKT');
面 百度
说明下：将wgs批量转到百度地图上，需要先转到百度经纬，再从百度经纬转到百度墨卡托，这些步骤被WGS2BDMKT封装了。
百度很麻烦很麻烦！！！
其他底图的偏移操作是，将wgs偏移到对应的底图，然后使用标准的墨卡托投影函数叠加到底图上即可。
百度底图操作是，将wgs偏移到对应的底图，但是底图是墨卡托坐标系，百度经纬度不能使用标准的墨卡托投影函数，二次加密的，所以要写个自定义函数，将百度经纬度自己转到百度墨卡托，才能叠加到百度底图。
这是其他转换库都没注意到，没涉及到的，坑爹的百度！！！
