矢量瓦片

GIS中的矢量与栅格数据

矢量切片(Vector tile)

在GIS中的数据分类有很多种方式，其中最常用的一种是根据数据组织结构方式的不同而分类成矢量数据和栅格数据的两种类型。

其中栅格数据以二维矩阵的形式来表示地理空间信息的数据结构，其中数据的最小存在单元是以像素的形式存在，可以理解为和图片的组织结构类似，以分辨率等特征作为精度的定义标准。

而矢量数据则是试图利用点、线、面等几何要素来表现这个世界，其数据结构紧凑精准，数据图形质量好，有利于地理信息检索与网络传输等。其中矢量数据的最小单元是以点的形式存在，点构成线，线组成面，面构造出体。所以，我个人看来矢量数据应该更贴近于信息的精准分析与计算，而栅格数据则偏重于信息的表达（主要受制于当前图像处理技术的瓶颈）。

栅格数据	矢量数据数据

为了带来更好更快的用户体验，目前许多主流WebGIS应用都采用了栅格切片技术，通过缓存切片的形式使得地图数据的浏览体验更顺畅。打开浏览器，F12出控制台，进入任意一家地图供应商提供的地图应用，你会发现大部分的地图数据都是以切片的形式请求获得的。栅格地图的切片应用是很广泛，可在我们的日常工作中遇到的需求往往要比这些功能需求较浅的商业性地图复杂，有的时候用户甚至会提出需要地图配色的编辑修改功能这样的需求，这是商业主流地图所达不到的，因为栅格切片在完成切图之后，所能控制的最小单位是一张图片，失去了对图片上地理信息的交互能力。

总结来看，栅格切片存在以下的几个缺点：

地图数据一次渲染，无法修改
无交互能力

那可否用WFS来替代呢？直接用WFS请求获取矢量数据，这样不就获得了交互能力吗？当然，如果在你的应用中矢量数据量不大的情况下，这样做也是可行的，但是一旦数据量大了起来，前端对于数据的请求和响应处理渲染会提高客户端的硬件门槛，而频繁的交互操作也会对服务器产生压力。

直接加载的矢量数据与对栅格地图进行切片这两种方式看起来好像有些互补，如果能将这二者结合起来的话应该会很美好：矢量+切片=矢量切片。

什么是矢量切片？

和栅格切片一样的思路，以金字塔的方式切割矢量数据，只不过切割的不是栅格图片，而是矢量数据的描述性文件，目前矢量切片主要有以下三种格式：GeoJSON，TopoJSON和MapbBox Vector Tile(MVT)。

栅格切图后文件存储形式：

矢量切图后文件存储形式：

切片中的数据结构：

从上面的两张图可以看出，其实思路是一致的，因此，矢量切图结合了矢量数据与栅格切图的优势互补：

前端缓存切片，提高地图使用的体验
粒度上来看，矢量切图继承了矢量数据的特性，以要素为单位进行管理，加强了细节上的把控能力
在保证体验的前提下，为用户提供地图数据样式动态修改的功能，加强了地图定制化的程度
数据的实时性

如何生成矢量切片？

矢量切片生成方式共有以下几种：

1）ArcGIS 系列产品：利用ArcGIS Pro生成矢量切片，然后发布在ArcGIS Online上；

2）Mapbox，目前已经提出了一套开放的矢量切片标准，并被多个开源团队所接受；

3）GeoServer，在2.11beta版中出现了对矢量切片的支持，主要依赖于开源插件geoserver-2.11-SNAPSHOT-vectortiles-plugin以及内嵌的GeoWebcahce完成切片工作。

==本文采取Mapbox的方式==

Mapbox的矢量瓦片

MVT简介

Mapbox的矢量瓦片（mapbox vector tile）是一种轻量级的数据格式，用于存储地理空间矢量数据，例如点、线和多边形。

矢量瓦片以 Google Protobufs（PBF）进行编码，这允许对结构化数据进行序列化。

Mapbox矢量瓦片使用的是.mvt的文件后缀。

Mapbox矢量瓦片标准

Mapbox矢量瓦片标准

Vector tiles standards

几何编码

Encoding geometry

几何编码将点、线和多边形编码为相对于网格左上角的x/y坐标对，如下图所示

MoveTo：起点

LineTo：画线

ClosePath：闭环

属性编码

Encoding attributes

属性编码是采用键值对的方式，压缩重复的属性值，将属性的键值对与矢量瓦片的键值对一一对应

有了矢量瓦片标准，对于在Web里面来说，利用WebGL的moveTo/lineTo之类的API，可以将矢量瓦片绘制出来。

使用PostGIS生成mvt

WebMecator投影瓦片计算规则

在生成mvt前，需要先了解一下瓦片的计算规则

参考链接：

https://zhuanlan.zhihu.com/p/548171000

https://www.cnblogs.com/beniao/archive/2010/04/18/1714544.html

墨卡托投影（Mercator Projection），又名“等角正轴圆柱投影”，荷兰地图学家墨卡托（Mercator）在1569年拟定，假设地球被围在一个中空的圆柱里，其赤道与圆柱相接触，然后再假想地球中心有一盏灯，把球面上的图形投影到圆柱体上，再把圆柱体展开，这就是一幅标准纬线为零度（即赤道）的“墨卡托投影”绘制出的世界地图。

由于赤道半径为6378137米，则赤道周长为2 * PI * r = 40075016.685578486162，因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。当纬度φ接近两极，即90°时，Y值趋向于无穷。因此通常把Y轴的取值范围也限定在[-20037508.3427892,20037508.3427892]之间。因此在墨卡托投影坐标系（米）下的坐标范围是：最小为(-20037508.3427892, -20037508.3427892 )到最大坐标为(20037508.3427892, 20037508.3427892)。

基础值

1
2
3
4


地球半径：EARTH_RADIUS = 6378137; 
X和Y的最大值： worldMercMax = π*EARTH_RADIUS = 20037508.3427892
X和Y的最小值： worldMercMin = -1 * worldMercMax
地球赤道周长： worldMercSize = worldMercMax - worldMercMin= 2 * π * EARTH_RADIUS = 40075016.68557848616

谷歌的瓦片地图规范已经实行多年，并深受认可，Bing、OSM、高德、天地图等都使用谷歌的瓦片地图规范，谷歌瓦片地图规范坐标原点为东经180°，北纬85.05°，==x轴向右，y轴向下==。

瓦片被切成256*256像素的图片，在不同的zoom等级，

t瓦片数量为： $2^{level}$

例如，zoom level为3时，横轴和纵轴方向瓦片数量均为8，瓦片对应坐标如下图所示。

生成mvt的ST函数

参考链接：

https://www.jianshu.com/p/e15ca286546c

https://blog.csdn.net/qq_35241223/article/details/106439268

主要函数：

ST_AsMvtGeom 用于将一个图层中位于参数box2d范围内的一个几何图形的所有坐标转换为MapBox Vector Tile坐标空间里的坐标。

ST_AsMVT 用于将基于MapBox Vector Tile坐标空间的几何图形转换为MapBox VectorTile二进制矢量切片

ST_TileEnvelope（3.0以上支持）根据行列号获取Envelope

辅助函数：

ST_Transform 坐标转换函数，用它可以做到支持任何坐标系的矢量瓦片

ST_Simplify 简化，用它来做线或者面的简化

ST_SimplifyPreserveTopology 与简化类似

ST_AsMvtGeom 应该将geom转到墨卡托坐标下，4326坐标下制作矢量瓦片有横向的压扁

ST_AsMVT

用于将基于MapBox Vector Tile坐标空间的几何图形转换为MapBox VectorTile二进制矢量切片

row —— 至少具有一个geometry列的行数据。
name —— 图层名字，默认为"default"。
extent —— 由MVT规范定义的屏幕空间（MVT坐标空间）中的矢量切片范围。
geom_name —— row参数的行数据中geometry列的列名，默认是第一个*geometry类型的列。
feature_id_name —— 行数据中要素ID列的列名。如果未指定或为NULL，则第一个有效数据类型（smallint, integer, bigint）的列将作为要素ID列，其他的列作为要素属性列。

ST_AsMVTGeom

用于将一个图层中位于参数box2d范围内的一个几何图形的所有坐标转换为MapBox Vector Tile坐标空间里的坐标。

geom —— 被转换的几何图形信息。
bounds—— 某个矢量切片的范围对应的空间参考坐标系中的几何矩形框（没有缓冲区）。
extent—— 是按规范定义的矢量切片坐标空间中的某个矢量切片的范围。如果为NULL，则默认为4096（边长为4096个单位的正方形）。
buffer—— 矢量坐标空间中缓冲区的距离，位于该缓冲区的几何图形部位根据clip_geom参数被裁剪或保留。如果为NULL，则默认为256。
clip_geom—— 用于选择位于缓冲区的几何图形部位是被裁剪还是原样保留。如果为NULL，则默认为true。

PostGIS生成MVT矢量切片的步骤

使用ST_AsMVTGeom函数将几何图形的所有坐标转换为MapBox Vector Tile坐标空间里的坐标，这样就将基于空间坐标系的几何图形转换成了基于MVT坐标空间的几何图形。
使用ST_AsMVT函数将基于MVT坐标空间的几何图形转换为MVT二进制矢量切片。

可以通过"||“操作符调用多次这个函数来同时创建多个图层的同一位置的矢量切片。

PostGIS ST函数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11


WITH 
	bounds AS ( SELECT 
		ST_Segmentize(ST_MakeEnvelope(13071343.332991, 3678761.297309, 13149614.849955, 3757032.814273, 3857),19567.879241) AS geom, 
		ST_Segmentize(ST_MakeEnvelope(13071343.332991, 3678761.297309, 13149614.849955, 3757032.814273, 3857),19567.879241)::box2d AS b2d
	), 
	mvtgeom AS ( SELECT ST_AsMVTGeom(ST_Transform(t.geom, 3857), bounds.b2d) AS geom, *  
	  FROM js_city_point_u_6492c1048429434b637f1be5 t, bounds 
	  WHERE ST_Intersects(t.geom, ST_Transform(bounds.geom, 4326))
	)
SELECT ST_AsMVT(mvtgeom.* , 'js_city_point_u_6492c1048429434b637f1be5' ) 
FROM mvtgeom

这段 PostGIS SQL 语句就是用于生成 MVT（Mapbox Vector Tiles）数据。

让我们逐步解释这段 SQL 语句的意思：

首先，使用 WITH 子句创建一个名为 bounds 的临时表，其中包含一个几何对象和一个边界框（box2d）。这个几何对象是通过使用 ST_MakeEnvelope 函数创建的，该函数根据提供的坐标范围**（tileToEnvelope函数计算得到）**和坐标系生成一个矩形。然后，使用 ST_Segmentize 函数对几何对象进行分段，以确保生成的矢量瓦片在显示时具有一定的细节和精度。
接下来，使用 mvtgeom 作为名称创建另一个临时表。在这个临时表中，使用 ST_Transform 函数将 t.geom（js_city_point_u_6492c1048429434b637f1be5 表中的几何列）从源坐标系（4326）转换为目标坐标系（3857）。然后，使用 ST_AsMVTGeom 函数将转换后的几何对象与 bounds.b2d 边界框进行相交（使用 ST_Intersects 函数），从而获得符合条件的几何对象。在结果中，还包含原始表 t 的所有列。
最后，使用 SELECT 语句从 mvtgeom 表中选择几何对象，并使用 ST_AsMVT 函数将这些几何对象转换为 MVT 格式的瓦片。'js_city_point_u_6492c1048429434b637f1be5' 是指定瓦片的图层名称。

简而言之，这个 SQL 查询的目的是生成包含符合条件的几何对象的 MVT 瓦片数据。它利用了 PostGIS 的空间函数和转换功能来处理几何对象的坐标系转换和分段，然后将结果转换为 MVT 格式的瓦片数据。

ST_MakeEnvelope需要的范围是如何确定的呢？

有了前面瓦片的计算规则的基础，就不难理解所求瓦片的x/y最大值最小值范围的计算方法了

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


HashMap<String, Double> tileToEnvelope(int zoom, int x, int y) {
    double worldMercMax = 20037508.3427892;
    double worldMercMin = -1 * worldMercMax;
    double worldMercSize = worldMercMax - worldMercMin;

    double worldTileSize = Math.pow(2, zoom);
    double tileMercSize = worldMercSize / worldTileSize;

    HashMap<String, Double> env = new HashMap<String, Double>();
    env.put("xmin", worldMercMin + tileMercSize * x);
    env.put("xmax", worldMercMin + tileMercSize * (x + 1));
    env.put("ymin", worldMercMax - tileMercSize * (y + 1));
    env.put("ymax", worldMercMax - tileMercSize * y);
    return env;
}

Java代码实现

接口返回的数据是==MVT二进制矢量切片==

Controller

1
2
3
4
5
6
7
8
9


@GetMapping(value = "/mvt/{tableName}/{zoom}/{x}/{y}.pbf")
public void getMvt(@PathVariable("tableName") String tableName,
                   @PathVariable("zoom") int zoom,
                   @PathVariable("x") int x,
                   @PathVariable("y") int y,
                   HttpServletResponse response) {

    pgService.getMvt(zoom, x, y, tableName, response);
}

Service

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32


public void getMvt(int zoom, int x, int y, String tableName, HttpServletResponse response) {
    try {
        String sql = getMvtSql(zoom, x, y, tableName);
        if (sql == null) {
            return;
        }
        log.info("DefaultPgSource: " + zoom + ", " + x + ", " + y + ":" + sql);

        byte[] mvtByte = mvtRepository.getMvtFromDefaultPg(sql);
        returnMvtByte(mvtByte, zoom, x, y, response);
    } catch (Exception e) {
        log.error(e.getMessage());
    }
}

public String getMvtSql(int zoom, int x, int y, String tableName) {
    if (!MvtUtils.tileIsValid(zoom, x, y)) {
        return null;
    }
    HashMap<String, Double> envelope = MvtUtils.tileToEnvelope(zoom, x, y);
    String sql = MvtUtils.envelopeToSQL(envelope, tableName);
    return sql;
}

public void returnMvtByte(byte[] mvtByte, int zoom, int x, int y, HttpServletResponse response) throws IOException {
    response.setHeader("Access-Control-Allow-Origin", "*");
    response.setHeader("Content-type", "application/vnd.mapbox-vector-tile");
    String mtvFileName = String.format("%d_%d_%d.mvt", zoom, x, y);
    response.setHeader("Content-Disposition", "attachment;filename=" + new String(mtvFileName.getBytes("UTF-8"), "iso-8859-1"));
    OutputStream os = response.getOutputStream();
    os.write(mvtByte);
}

Repository

1
2
3
4
5
6
7
8
9


public byte[] getMvtFromDefaultPg(String sql) {
    try {
        byte[] reByte = jdbcTemplate.queryForObject(sql, (rs, rowNum) -> rs.getBytes("st_asmvt"));
        return reByte;
    } catch (Exception e) {
        log.error("默认数据库瓦片获取失败" + e.getMessage());
        return null;
    }
}

MVT工具类

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58


public class MvtUtils {

    public static Boolean tileIsValid(int zoom, int x, int y) {
        Double size = Math.pow(2, zoom);
        if (x >= size || y >= size) {
            return false;
        }
        if (x < 0 || y < 0) {
            return false;
        }
        return true;
    }

    public static HashMap<String, Double> tileToEnvelope(int zoom, int x, int y) {
        double worldMercMax = 20037508.3427892;
        double worldMercMin = -1 * worldMercMax;
        double worldMercSize = worldMercMax - worldMercMin;

        double worldTileSize = Math.pow(2, zoom);
        double tileMercSize = worldMercSize / worldTileSize;

        HashMap<String, Double> env = new HashMap<String, Double>();
        env.put("xmin", worldMercMin + tileMercSize * x);
        env.put("xmax", worldMercMin + tileMercSize * (x + 1));
        env.put("ymin", worldMercMax - tileMercSize * (y + 1));
        env.put("ymax", worldMercMax - tileMercSize * y);
        return env;
    }

    public static String envelopeToBoundsSQL(HashMap<String, Double> env) {
        int DENSIFY_FACTOR = 4;
        env.put("segSize", (env.get("xmax") - env.get("xmin")) / DENSIFY_FACTOR);
        String sqlTemp = String.format("ST_Segmentize(ST_MakeEnvelope(%f, %f, %f, %f, 3857),%f)", env.get("xmin"), env.get("ymin"), env.get("xmax"), env.get("ymax"), env.get("segSize"));
        return sqlTemp;
    }

    public static String envelopeToSQL(HashMap<String, Double> env, String tableName) {
        // lines_pg   gis_osm_transport_free_1
        HashMap<String, String> table = new HashMap<String, String>();
        table.put("table", tableName);
        table.put("srid", "4326");
        table.put("geomColumn", "geom");
        table.put("attrColumns", " * ");
        table.put("env", envelopeToBoundsSQL(env));

        String mvtsql = MessageFormat.format("WITH" +
                " bounds AS ( SELECT {0} AS geom, {0}::box2d AS b2d)," +
                " mvtgeom AS (" +
                " SELECT ST_AsMVTGeom(ST_Transform(t.{1}, 3857), bounds.b2d) AS geom, {2}" +
                " FROM {3} t, bounds" +
                " WHERE ST_Intersects(t.{1}, ST_Transform(bounds.geom, {4}))" +
                " )" +
                " SELECT ST_AsMVT(mvtgeom.* , ''{3}'' ) FROM mvtgeom",
            table.get("env"), table.get("geomColumn"), table.get("attrColumns"), table.get("table"), table.get("srid"));

        return mvtsql;
    }
}

矢量瓦片#

GIS中的矢量与栅格数据#

什么是矢量切片？#

如何生成矢量切片？#

Mapbox的矢量瓦片#

MVT简介#

Mapbox矢量瓦片标准#

几何编码#

属性编码#

使用PostGIS生成mvt#

WebMecator投影瓦片计算规则#

生成mvt的ST函数#

ST_AsMVT#

ST_AsMVTGeom#

PostGIS生成MVT矢量切片的步骤#

PostGIS ST函数#

Java代码实现#

Controller#

Service#

Repository#

MVT工具类#