java+postgis实现4490坐标系动态矢量瓦片并实现后端聚合

前言

什么是矢量瓦片：现阶段，电子地图瓦片主要使用两种方式，一种是传统的栅格瓦片，另外一种是新出的矢量瓦片（Vector Tiles），前者是采用四叉树金字塔模型的分级方式，将地图切割成无数大小相等的矩形栅格图片，由这些矩形栅格图片按照一定规则拼接成不同层级的地图显示。后者是将矢量数据用多层次模型分割成矢量要素描述文件存储在服务器端，再到客户端根据指定样式进行渲染绘制地图，在单个矢量瓦片上存储着投影于一个矩形区域内的几何信息和属性信息。当客户端通过分布式网络获取矢量瓦片、地图标注字体、图标、样式文件等数据后，最终在客户端进行渲染输出地图。

简而言之：矢量瓦片存储的实际还是矢量数据优点数据不会丢失精度，同时在一定程度上也解决了矢量数据加载和渲染慢问题，是目前最流行的一种地图数据解决方案。

这里介绍下为什么要4490坐标系呢？

其实地理数据的坐标系是多种多样的，比如常见的全球常用的通用坐标系WGS84和百度的BD09和高德的GJC-02,随着国家技术发展日新月异，国产化潮流来袭，EPSG:4490坐标系成为了目前国内项目上逐渐的地图坐标系标准,CGCS2000(4490) : 国家大地坐标系。天地图采用的就是这个。可以提供高精度三维坐标，它是以地球质量中心为原点的地心坐标系，可以满足航天、海洋、气象、水利、建设、规划、地质调查、国土资源管理等领域的多种需求。所以后面地理信息系统使用4490坐标系的要求只会越来越严格。

一、postgis是什么？

postgis是基于postgresql数据库的一款空间插件，PostgreSQL目前使用量也是越来越多逐渐赶超主流的mysql等数据库，postgis在这块提供里非常丰富的空间数据处理函数和扩展，是目前gis开发人员使用最多的数据库和插件。

二、本文主要使用的postgis函数

ST_AsMVTGeom：将几何图形转换成MAPBOX矢量瓦片的坐标空间

参数介绍：
geom:指定需要变换的几何字段
bounds:不带缓冲区的几何边界
extent:按照规范定义的平铺坐标空间中的平铺范围。默认4096
buffer:平铺坐标空间中任意修剪几何图形的缓冲区距离,默认256
clip_geom:用于控制几何图形是否被裁剪活编码，默认true

ST_AsMVT：返回一个MapBox矢量瓦片的一组行（返回Byte数组）

参数介绍：
row:一个几何列的行数据
name:图层名称
extent:按照规范定义的平铺坐标空间中的平铺范围。默认4096 （非必填）
geom_name:指定行数据列里面的空间字段名称，默认是找到的第一个空间字段 （非必填）

ST_MakeEnvelope：从X和Y的最小值和最大值创建矩形。输入值必须在SRID指定的空间参考系中。如果未指定SRID，则使用未知空间参考系(SRID 0)。

参数介绍：
geometry ST_MakeEnvelope(float xmin, float ymin, float xmax, float ymax, integer srid=unknown);
这里就是数据的四至范围，最后一个参数为坐标系的srid，这里是4490

st_centroid：计算几何的几何中心，或等效地将几何的质心计算为POINT（后端聚合使用）

st_collect：从其他Geometry对象的collection返回一个具体的Geometry值（后端聚合使用）

width_bucket：将对应的经度、纬度范围值划分为多少份，参数就是地图四至范围和划分数量（后端聚合使用）

st_intersects：判断两个几个对象是否存在相交，返回布尔值。（数据过滤常用）

st_x：获取几何的x坐标值
st_y：获取几何的y坐标值

三、实现4490坐标系动态切片

1.核心代码

代码如下（示例）：

        String tile = TileUtils4326.xyz2prj4326wkt(z, x, y);
        String sql = "select ST_AsMVT ( mvtgeom.*, '" + sourceName + "' )as data from ( " +
                "   SELECT " +
                "      ST_AsMVTGeom ( the_geom, ST_MakeEnvelope ( "+TileUtils4326.xyz2prj4326BBox(z, x, y)+", 4490 ),4096,0,true ) AS the_geom, " +
                "      * " +
                "   FROM " +
                "      \""+sourceName+"\" " +
                "   where st_intersects(the_geom,st_geomfromtext(?,4326))"+
                "    )mvtgeom";
        Map<String, Object> result = jdbcTemplate.queryForMap(sql, tile);
        ByteArrayOutputStream out = new ByteArrayOutputStream();
        GZIPOutputStream gzip;
        try {
            gzip = new GZIPOutputStream(out);
            gzip.write((byte[]) result.get("data"));
            gzip.close();
        } catch (IOException e) {
            logger.error("gzip compress error.", e);
        }
        return out.toByteArray();

2.xyz行列号转经纬度工具类

代码如下（示例）：

public class TileUtils4326 {
    public static String xyz2prj4326wkt(int z, int x, int y) {
        StringBuilder sb = new StringBuilder("POLYGON ((");
        double n = Math.pow(2, z);
        double lon_min = (x / n) * 360.0 - 180.0;
        double lat_min = 90.0 - (((y + 1) / n) * 360.0);
        double lon_max = ((x + 1) / n) * 360.0 - 180.0;
        double lat_max = 90.0 - ((y / n) * 360.0);
        sb.append(lon_min +" "+lat_max+", ");
        sb.append(lon_max +" "+lat_max+", ");
        sb.append(lon_max +" "+lat_min+", ");
        sb.append(lon_min +" "+lat_min+", ");
        sb.append(lon_min +" "+lat_max+")) ");
        return sb.toString();
    }

    public static String xyz2prj4326BBox(int z, int x, int y) {
        String bbox = "";
        double n = Math.pow(2, z);
        double lon_min = (x / n) * 360.0 - 180.0;
        double lat_min = 90.0 - (((y + 1) / n) * 360.0);
        double lon_max = ((x + 1) / n) * 360.0 - 180.0;
        double lat_max = 90.0 - ((y / n) * 360.0);
        bbox = lon_min + ","+lat_min+","+lon_max+","+lat_max;
        return bbox;
    }
}

四、实现后端聚合4490坐标系动态切片

1.核心代码

bbox为数据的四至范围

String tile = TileUtils4326.xyz2prj4326wkt(z, x, y);
String sql = "select ST_AsMVT ( mvtgeom.*, '" + sourceName + "' )as data from (SELECT " +
        "width_bucket(st_x(the_geom), " + bbox[0] + " ," + bbox[2] + " ," + segmentedCount + ") grid_x," +
        "width_bucket(st_y(the_geom), " + bbox[1] + " ," + bbox[3] + " ," + segmentedCount + ") grid_y," +
        "count(*)," +
        "ST_AsMVTGeom (st_centroid(st_collect(the_geom)), ST_MakeEnvelope ( " + TileUtils4326.xyz2prj4326BBox(z, x, y) + ", 4490 ),4096,0,true ) as geom " +
        "from \"" + sourceName + "\" " +
        "where " +
        "st_x(the_geom) between " + bbox[0] + " and " + bbox[2] + " " +
        "and " +
        "st_y(the_geom) between " + bbox[1] + "  and " + bbox[3] + " " +
        "and " +
        "st_intersects(the_geom,st_geomfromtext(?,4326)) " +
        "GROUP BY grid_x,grid_y )mvtgeom";
Map<String, Object> result = jdbcTemplate.queryForMap(sql, tile);
ByteArrayOutputStream out = new ByteArrayOutputStream();
GZIPOutputStream gzip;
try {
     gzip = new GZIPOutputStream(out);
     gzip.write((byte[]) result.get("data"));
     gzip.close();
} catch (IOException e) {
     logger.error("gzip compress error.", e);
}
return out.toByteArray();

2.效果图