使用PostgreSQL+PostGIS实现地图引擎（经历ArcGIS Server/GeoServer/MapServer后）_pg sql搭建地图服务

地图服务器方案探索

研究GIS服务器已经很长一段时间，最开始使用ArcGIS Server,后来由于费用问题公司基于nginx和lua自己做了一套矢量瓦片的地图服务器，再后来开始研究使用Geoserver,再后来由于性能问题改研究Mapserver，最近由于开发的灵活性和业务需求开始研究PostgreSQL+PostGIS的方案。

PostGIS官方方案

WITH mvtgeom AS
(
  SELECT ST_AsMVTGeom(geom, ST_TileEnvelope(12,513,412)) AS geom, name, description
  FROM points_of_interest
  WHERE ST_Intersects(geom, ST_TileEnvelope(12,513,412)
)
SELECT ST_AsMVT(mvtgeom.*)
FROM mvtgeom;
1
2
3
4
5
6
7
8

存在问题

ST_TileEnvelope只能生成3857坐标的边界，类似于：
POLYGON((11858134.820049 3561354.02186289,11858134.820049 3580921.90110389,11877702.69929 3580921.90110389,11877702.69929 3561354.02186289,11858134.820049 3561354.02186289))

我们数据库里的数据都是4490坐标的边界，即类似于：
POLYGON((-180 -90,-180 90,180 90,180 -90,-180 -90))

如果直接使用ST_TileEnvelope生成的边界来计算就必须要求数据库里的数据转换为3857坐标系的样式。

解决办法有五个思路，一是修改源数据为3857的数据后再导入，二是4490坐标系的数据设置坐标系为3857，三是在计算过程中用st_transform强制投影转换，四是计算过程中采用::geography强制转换，五是放弃ST_TileEnvelope函数采用ST_MakeEnvelope并从代码中计算好4490下的瓦片边界后使用。

方案一：
不可取。由于国家要求和实际计算的需要，我们只能在库里存4490坐标系的数据。

方案二：

不可取。
update "DS_YJJBNT_5" set geom=st_setsrid(geom,4490)这一行语句只是修改了表里面的坐标系描述，即只改了元数据，并不会去动表中的任何一条数据，所以设置坐标系的方式是无法使用的。

SELECT st_srid(geom) FROM "DS_YJJBNT_5" limit 1;查询数据表的坐标系。

方案三：
不可取。
官方文档中提示说使用st_transform需要有proj4这个库编译进PostgreSQL库，但是重新编译又是一项大工程。
UPDATE “DS_YJJBNT_5” SET geom=st_transform(st_geomfromtext((st_astext(geom)),4490),3857);

方案四：
不可取。
还没试成功过，具体原理未知。

方案五：
可以成功。

方案五的实现

通过xyz生成任意瓦片矩形边界

python版

func tile2lon( x int,  z int)(a float64) {
	return float64(x) /math.Pow(2, float64(z)) * 360.0 - 180;
 }

 func tile2lat( y int,  z int)(a float64) {
   n := math.Pi - (2.0 * math.Pi * float64(y)) / math.Pow(2, float64(z));
   return math.Atan(math.Sinh(n))*180/math.Pi;
 }

ymax :=FloatToString(tile2lat(int(xyz.y), int(xyz.z)));
ymin := FloatToString(tile2lat(int(xyz.y+1), int(xyz.z)));
xmin := FloatToString(tile2lon(int(xyz.x), int(xyz.z)));
xmax := FloatToString(tile2lon(int(xyz.x+1), int(xyz.z)));
1
2
3
4
5
6
7
8
9
10
11
12
13

node.js版

let xmin = x/Math.pow(2,z)*360.0-180;

let n = Math.PI - (2.0 * Math.PI * y) / Math.pow(2,z);
let ymax = Math.atan((Math.exp(n)-Math.exp(-n))/2)*180/Math.PI;
1
2
3
4

Sinh(n)与(Math.exp(n)-Math.exp(-n))/2等价

原版SQL

WITH mvtgeom AS (SELECT ST_AsMVTGeom(geom, ST_MakeEnvelope(106.5234375, 30.44867367928757, 106.69921875, 30.600093873550065, 4490)) AS geom, bsm, objectid 
FROM "DS_YJJBNT_5" WHERE ST_Intersects(geom, ST_MakeEnvelope(106.5234375, 30.44867367928757, 106.69921875, 30.600093873550065, 4490))) 
SELECT ST_AsMVT(mvtgeom.*,'points') FROM mvtgeom;
1
2
3

解决了地图缩放到任何层级切出来图片大小都会盖全图的问题

（设置extent=4096，buffer为0，裁剪为true）

WITH mvtgeom AS (SELECT ST_AsMVTGeom(geom, ST_MakeEnvelope(106.5234375, 30.44867367928757, 106.69921875, 30.600093873550065, 4490),4096, 0, true) AS geom, bsm, objectid 
FROM "DS_YJJBNT_5" WHERE ST_Intersects(geom, ST_MakeEnvelope(106.5234375, 30.44867367928757, 106.69921875, 30.600093873550065, 4490))) 
SELECT ST_AsMVT(mvtgeom.*,'points') FROM mvtgeom;
1
2
3

解决了放大到最大层级时相对geoserver部分数据丢失且偏移的问题

(去除where相交的条件）

WITH mvtgeom AS (SELECT ST_AsMVTGeom(geom, ST_MakeEnvelope(106.5234375, 30.44867367928757, 106.69921875, 30.600093873550065, 4490),4096, 0, true) AS geom, bsm, objectid 
FROM "DS_YJJBNT_5" ) 
SELECT ST_AsMVT(mvtgeom.*,'points') FROM mvtgeom;
1
2
3

重组美化SQL后

select ST_AsMVT ( mvtgeom.*, '${tableName}' )as data from (
   SELECT
      ST_AsMVTGeom ( geom, ST_MakeEnvelope ( #{xmin},#{ymin},#{xmax},#{ymax}, 4490 ),4096,0,true ) AS geom,
      bsm,
      objectid
   FROM
      "${tableName}"
    )mvtgeom
1
2
3
4
5
6
7
8

参考资料：

https://www.cnblogs.com/polong/p/9831106.html

http://postgis.net/docs/manual-3.0/ST_TileEnvelope.html

http://postgis.net/docs/manual-3.0/ST_AsMVT.html

http://postgis.net/docs/manual-3.0/ST_AsMVTGeom.html

我自己产品的官网：http://www.newgis.top