地图服务器方案探索
研究GIS服务器已经很长一段时间，最开始使用ArcGIS Server,后来由于费用问题公司基于nginx和lua自己做了一套矢量瓦片的地图服务器，再后来开始研究使用Geoserver,再后来由于性能问题改研究Mapserver，最近由于开发的灵活性和业务需求开始研究PostgreSQL+PostGIS的方案。
2021.2.1更新
最近实际使用中发现了些问题，在重庆范围的数据，叠加在3857的天地图上时，如果放大到14级以上看不出来偏移，便是缩小到能看全中国，能看全球的时候地图就偏移到了西伯利亚去了。
所以有条件还是应该选用“方案三”在地图生成过程中使用st_transform函数进行坐标转换。
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;
存在问题
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)));
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;
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;
解决了地图缩放到任何层级切出来图片大小都会盖全图的问题
（设置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;
解决了放大到最大层级时相对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;
重组美化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
GIS引擎开发
最近自己用Mapserver包装了一个GIS引擎作为Geoserver的替代方案，适合小微企业和个人用户，下载地址：
解压密码：2234，后续我将按计划进行完善，力争做到轻量、易用、稳定、高性能、开源。
同时我将按低价提供服务，并低价出售源代码，以保证日常开销。
参考资料：
https://www.cnblogs.com/polong/p/9831106.html
http://postgis.net/docs/manual-3.0/ST_TileEnvelope.html
