蜂巢的艺术与技术价值 - PostgreSQL PostGIS's hex-grid

背景

人们为了更好的描述一个东西，有一种将大化小的思路，比如时钟被分为了12个区域，每个区域表示一个小时，然后每个小的区域又被划分为更小的区域表示分钟。

在GIS系统中，也有类似的思想，比如将地图划分成网格。通过编码来简化地理位置的判断（比如相交，包含，距离计算等），但是请注意使用网格带来的问题，比如精度的问题，网格的大小决定了精度，又比如相对坐标的问题，可能无法描述清楚边界的归属。

PS：

1. 在PostGIS中虽然也支持网格对象的描述方式，但是并不是使用这种方法来进行几何运算（比如相交，包含，距离计算等），所以不存在类似的精度问题，个人建议没有强需求的话，不必做这样的网格转换。

PostgreSQL GIS索引的原理请参考

《从难缠的模糊查询聊开 - PostgreSQL独门绝招之一 GIN , GiST , SP-GiST , RUM 索引原理与技术背景》

2. 如果是多种精度地图的切换（比如多个图层，每个图层代表一种地图精度），建议使用辐射的方式逐渐展开更精细的图层，以点为中心，逐渐辐射。（很多专业的地图软件是这样做的）

回到主题，还记得最强大脑的蜂巢迷宫吗？

还有勤劳的蜜蜂兄弟

我们接下来要做的就是如何将几何图形转换为网格对象。

回忆一下六边形的几何特性

首先要了解一下六边形的几何特性，提供转换的计算基础。

六边形可以切分为6个等边三角形

所以它的边长关系如下

面积计算

将几何图形（sharp）转换为六边形网格

比如要将澳大利亚的版图，转换为六边形网格，

有两种方法，一种使用geotools的JAVA 类（在程序中转换），另一种是使用PostGIS插件的UDF（在数据库中转换）。

当然，如果PostgreSQL安装了pljava插件的话，那么也可以在PostgreSQL中调用geotools提供的java类进行转换。

下面是例子

1 geotools Vector grids class

http://docs.geotools.org/latest/userguide/extension/grid.html

使用geotools vector grids class生成网格，返回 SimpleFeatureSource 类型。

geotools Vector grids class支持将几何图形转换为 polygon网格或者 line网格。

1 Polygon grids

举几个例子

1. 将澳大利亚地图转换为10度边长的正方形网格

输入澳大利亚的经纬度范围，转换

    ReferencedEnvelope gridBounds = new ReferencedEnvelope(
            110.0, 150.0, -45.0, -5.0, DefaultGeographicCRS.WGS84);  

    SimpleFeatureSource grid = Grids.createSquareGrid(gridBounds, 10.0);

2. 将澳大利亚地图转换为最大20度边长的扇形网格

    ReferencedEnvelope gridBounds = new ReferencedEnvelope(
            110, 160, -45, -8, DefaultGeographicCRS.WGS84);  

    double squareWidth = 20.0;  

    // max distance between vertices
    double vertexSpacing = squareWidth / 20;  

    SimpleFeatureSource grid = Grids.createSquareGrid(gridBounds, squareWidth, vertexSpacing);

3. 创建纵横宽100，变长为5.0的六边形网格

    ReferencedEnvelope gridBounds = new ReferencedEnvelope(0, 100, 0, 100, null);  

    // length of each hexagon edge
    double sideLen = 5.0;
    SimpleFeatureSource grid = Grids.createHexagonalGrid(gridBounds, sideLen);

4. 导入图形、将澳大利亚地图转换为边长1度的六边形网格

自定义图形边界判断的类

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import java.io.IOException;
import java.util.Map;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.factory.CommonFactoryFinder;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.filter.Filter;
import org.opengis.filter.FilterFactory2;  

public class IntersectionBuilder extends GridFeatureBuilder {
    final FilterFactory2 ff2 = CommonFactoryFinder.getFilterFactory2();
    final GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();  

    final SimpleFeatureSource source;
    int id = 0;  

    public IntersectionBuilder(SimpleFeatureType type, SimpleFeatureSource source) {
        super(type);
        this.source = source;
    }  

    public void setAttributes(GridElement el, Map<String, Object> attributes) {
        attributes.put("id", ++id);
    }  

    @Override
    public boolean getCreateFeature(GridElement el) {
        Coordinate c = ((PolygonElement) el).getCenter();
        Geometry p = gf.createPoint(c);
        Filter filter = ff2.intersects(ff2.property("the_geom"), ff2.literal(p));
        boolean result = false;  

        try {
            result = !source.getFeatures(filter).isEmpty();
        } catch (IOException ex) {
            throw new IllegalStateException(ex);
        }  

        return result;
    }
}

导入地图，在createHexagonalGrid中使用边界判断的类，生成六边形网格

    // Load the outline of Australia from a shapefile
    URL url = getClass().getResource("oz.shp");
    FileDataStore dataStore = FileDataStoreFinder.getDataStore(url);
    SimpleFeatureSource ozMapSource = dataStore.getFeatureSource();  

    // Set the grid size (1 degree) and create a bounding envelope
    // that is neatly aligned with the grid size
    double sideLen = 1.0;
    ReferencedEnvelope gridBounds =
            Envelopes.expandToInclude(ozMapSource.getBounds(), sideLen);  

    // Create a feature type
    SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
    tb.setName("grid");
    tb.add(GridFeatureBuilder.DEFAULT_GEOMETRY_ATTRIBUTE_NAME,
            Polygon.class, gridBounds.getCoordinateReferenceSystem());
    tb.add("id", Integer.class);
    SimpleFeatureType TYPE = tb.buildFeatureType();  

    // Build the grid the custom feature builder class
    GridFeatureBuilder builder = new IntersectionBuilder(TYPE, ozMapSource);
    SimpleFeatureSource grid = Grids.createHexagonalGrid(gridBounds, sideLen, -1, builder);

5. 生成的六边形网格的摆放参数，横(flat)的还是竖(angled)的？

默认为flat网格，如果要生成angled网格，如下

    ReferencedEnvelope gridBounds = new ReferencedEnvelope(0, 100, 0, 100, null);
    double sideLen = 5.0;
    GridFeatureBuilder builder = new DefaultGridFeatureBuilder();
    SimpleFeatureSource grid = Hexagons.createGrid(
            gridBounds, sideLen, HexagonOrientation.ANGLED, builder);

2 line grids

转换为line 网格

例子

    ReferencedEnvelope gridBounds = new ReferencedEnvelope(
            110.0, 150.0, -45.0, -5.0, DefaultGeographicCRS.WGS84);  

    /*
     * Line definitions:
     * major lines at 10 degree spacing are indicated by level = 2
     * minor lines at 2 degree spacing are indicated by level = 1
     * (level values are arbitrary; only rank order matters)
     */
    List<OrthoLineDef> lineDefs = Arrays.asList(
            // vertical (longitude) lines
            new OrthoLineDef(LineOrientation.VERTICAL, 2, 10.0),
            new OrthoLineDef(LineOrientation.VERTICAL, 1, 2.0),  

            // horizontal (latitude) lines
            new OrthoLineDef(LineOrientation.HORIZONTAL, 2, 10.0),
            new OrthoLineDef(LineOrientation.HORIZONTAL, 1, 2.0));  

    // Specify vertex spacing to get "densified" polygons
    double vertexSpacing = 0.1;
    SimpleFeatureSource grid = Lines.createOrthoLines(gridBounds, lineDefs, vertexSpacing);

2 PostGIS UDF hex-grid

https://github.com/minus34/postgis-scripts/tree/master/hex-grid

PostGIS不需要多介绍了，几十年的老牌GIS插件，在军方、科研、民用等各个领域有着非常广泛对应用。

如果你使用了PostGIS插件的话，在里面存储了不管是geometry, polygon还是其他的地图类型，都可以转换为六边形网格。

转换时使用这些定义好的UDF即可。

UDF使用方法

See the 2 sample usage scripts to see how to create a national hex grid, using the function.

UDF输入参数

Parameter	Description
areakm2	Area of each hexagon in square km. Note: output hexagon sizes can be off slightly due to coordinate rounding in the calcs.
xmin,ymin	Minimum coordinates of the grid extents (i.e. bottom, left).
xmax,ymax	Maximum coordinates of the grid extents (i.e. top, right).
inputsrid	The coordinate system (SRID) of the min/max coordinates.
workingsrid	The SRID used to process the hexagons:
-	SRID must be a projected coord sys (i.e. in metres) as the calcs require ints. Degrees are out.
-	Should be an equal area SRID - i.e. Albers or Lambert Azimuthal (e.g. Australia = 3577, US = 2163).
-	Using a Mercator projection will NOT return hexagons of equal area (don't try it in Greenland).
ouputsrid	The SRID of the output hexagons.