使用postgis发布mvt矢量切片

若想发布自己矢量切片，除了使用mapbox推荐的工具与geoserver外还能够直接使用postgis+后端程序进行动态发布。

使用postgis动态发布矢量切片，必须使用postgis 2.4.0以上版本。需要使用以下两个函数：

ST_AsMVT
Synopsis
bytea ST_AsMVT(anyelement set row);

bytea ST_AsMVT(anyelement row, text name);

bytea ST_AsMVT(anyelement row, text name, integer extent);

bytea ST_AsMVT(anyelement row, text name, integer extent, text geom_name);

Description
Return a Mapbox Vector Tile representation of a set of rows corresponding to a Layer. Multiple calls can be concatenated to a tile with multiple Layers. Geometry is assumed to be in tile coordinate space and valid as per specification. Typically ST_AsMVTGeom can be used to transform geometry into tile coordinate space. Other row data will be encoded as attributes.

The Mapbox Vector Tile format can store features with a different set of attributes per feature. To make use of this feature supply a JSONB column in the row data containing Json objects one level deep. The keys and values in the object will be parsed into feature attributes.
row row data with at least a geometry column.

name is the name of the Layer. If NULL it will use the string "default".

extent is the tile extent in screen space as defined by the specification. If NULL it will default to 4096.

geom_name is the name of the geometry column in the row data. If NULL it will default to the first found geometry column.

Enhanced: 2.5.0 - added support parallel query.

Availability: 2.4.0


ST_AsMVTGeom

Synopsis
geometry ST_AsMVTGeom(geometry geom, box2d bounds, integer extent=4096, integer buffer=256, boolean clip_geom=true);

Description
Transform a geometry into the coordinate space of a Mapbox Vector Tile of a set of rows corresponding to a Layer. Makes best effort to keep and even correct validity and might collapse geometry into a lower dimension in the process.

geom is the geometry to transform.

bounds is the geometric bounds of the tile contents without buffer.

extent is the tile extent in tile coordinate space as defined by the specification. If NULL it will default to 4096.

buffer is the buffer distance in tile coordinate space to optionally clip geometries. If NULL it will default to 256.

clip_geom is a boolean to control if geometries should be clipped or encoded as is. If NULL it will default to true.

Availability: 2.4.0

发布矢量切片的逻辑如下：

graph LR
xyz-->瓦片空间范围
瓦片空间范围-->检索空间表
检索空间表-->转化成mvt格式文件

在这个流程中，后面三步都postgis能够帮助实现。因此关键就是将x y z转化成该瓦片的空间范围。这一步可以在后端实现、也可以在数据库实现。在此处通过数据库函数实现这个功能。

CREATE OR REPLACE FUNCTION BBox(x integer, y integer, zoom integer)
    RETURNS geometry AS
$BODY$
DECLARE
    max numeric := 6378137 * pi();
    res numeric := max * 2 / 2^zoom;
    bbox geometry;
BEGIN
    return st_transform(ST_MakeEnvelope(
        -max + (x * res),
        max - (y * res),
        -max + (x * res) + res,
        max - (y * res) - res,
        3857),4326);
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE;

SELECT ST_AsMVT(q, 'admin', 4096, 'geom')
FROM (
  SELECT id, name, admin_level,
    ST_AsMvtGeom(
      geometry,
      BBox(16597, 11273, 15),
      4096,
      256,
      true
    ) AS geom
  FROM import.osm_admin
  WHERE geometry && BBox(16597, 11273, 15)
  AND ST_Intersects(geometry, BBox(16597, 11273, 15))
) AS q;

这两个脚本就是核心代码。
下面使用Python写后台代码：

import psycopg2
#该函数的作用就是根据表名，动态生成sql语句
def createSql(table_name,x,y,z):
    sql1="SELECT (select f_geometry_column from geometry_columns where f_table_name= '"+table_name+"' ),replace(replace(string_agg(a.attname,','),(select f_geometry_column from geometry_columns where f_table_name= '"+table_name+"' )||',',''),(select f_geometry_column from geometry_columns where f_table_name='minemap_waterface')||'','') as name FROM pg_class as c,pg_attribute as a  where c.relname = '"+table_name+"' and a.attrelid = c.oid and a.attnum>0";
    conn = psycopg2.connect(database="postgis_db", user="postgres", password="postgres", host="192.168.25.95", port="5432");
    cur = conn.cursor();
    cur.execute(sql1);
    keyData = cur.fetchall();
    geom_column=keyData[0][0];
    other_fields=keyData[0][1];
    cur.close();
    conn.close();
    sql2="SELECT ST_AsMVT(q, '"+table_name+"', 4096, 'geometry') from (SELECT "+other_fields+", ST_AsMvtGeom("+geom_column+",BBox("+str(x)+", "+str(y)+" ,"+str(z)+") , 4096,256, true) as geometry  FROM "+table_name+" WHERE "+geom_column+" && BBox("+str(x)+", "+str(y)+" ,"+str(z)+")  AND ST_Intersects("+geom_column+", BBox("+str(x)+", "+str(y)+" ,"+str(z)+") )) q;"
    return sql2;

#运行sql,获得mvt的16进制文件
def getPbfFile(table_name,x,y,z):
    sql=createSql(table_name,x,y,z);
    conn = psycopg2.connect(database="postgis_db", user="postgres", password="postgres", host="192.168.25.95",
                            port="5432");
    cur = conn.cursor();

    cur.execute(sql);
    print(sql)
    keyData = cur.fetchall();
    cur.close();
    conn.close();
    with open(str(z)+"_"+str(x)+"_"+str(y)+'.pbf', 'wb') as f:
        f.write(keyData[0][0])
    return keyData[0][0];

tile=getPbfFile("minemap_waterface_copy",54605,26956,16);
print(tile)

G8Y%U%XA1SRX3(7O`5MIVHP.png

从文件结果来看，结果是对的。但是这个代码还是有问题的，

1、这样发布图层还是不太灵活。
2、取一次文件需要两次运行sql语句效率还是比较低的。
3、不太方便将多个表格的数据融合在一起，发布成一个图层。

使用postgis发布mvt矢量切片

推荐阅读更多精彩内容