GIS综合应用：高速公路结冰预警

前言

结冰对高速公路的安全行车具有很大的影响，在寒冷季节，温度低于0度的道路容易出现结冰现象，导致行车事故的发生。为了降低事故风险，对高速公路的结冰情况进行实时监测和预警分析具有重要意义。本项目使用了一种基于 GIS 和插值方法的全速公路结冰情况预测方法。

本项目仅考虑了温度对结冰的影响，实际上结冰还受到其他因素（如路面湿度、风速等）的影响，如果读者想获取更准确的预测结果，可以结合其他气象要素进行综合分析，进一步完善模型。

本项目使用全开源的技术栈与中间件，符合国产化要求。具体技术栈如下图：

空间处理	geotools
可视化	openlayer.js
空间数据库	postgis
地图服务器	geoserver

摘要

本项目基于 GIS 和插值方法的全国高速公路结冰情况预测。首先，收集气象站点的温度数据和高速公路的 GIS 数据。然后，利用插值方法计算温度等值面，并提取温度低于0度的区域。最后，通过相交分析得出潜在结冰的高速公路区域，并进行可视化和预警发布。可为高速公路结冰情况的实时监测和预警提供了一种有效的手段。

流程

数据收集与整理

数据准备

1）周期性收集了气象站点的预报数据，包括地理位置（经纬度）、气温、预警时间等。

2）收集了高速公路的 GIS 数据，包括公路的地理位置。

3）收集预警区域的GIS数据。

数据皆来源于互联网。

图层名称
自动站点最低温图层	station_tmin

构建图层

在 PostGIS 数据库中创建一个新的空间表，用于存储自动站点数据。表结构应包含以下字段：

• id：唯一标识符，例如整数或 UUID。
• station_name：站点名称。
• temperature：温度。
• geom：点几何对象，用于存储站点的空间位置信息。该字段应设置为“GEOMETRY(Point, )”类型，其中 是空间参考系统的代码。
• report_time：预测时间。

数据导入

使用 GeoTools 对 PostGIS 进行操作，将气象站点数据导入到新建的空间表中。

数据清洗

对气象站点数据进行预处理，清洗异常和缺失值。

插值方法与温度等值面计算

温度插值

利用克里金插值等插值方法，将气象站点的温度数据扩展到整个预警区域。

克里金插值法是一种广泛应用于空间插值的方法，它可以通过已知点的气象数据推断未知位置的气象数据，并且能够在保持较高插值精度的同时，有效地处理不规则和空间分布不均匀的气象数据。

构建栅格

读取预警区域的空间信息，构建1000*1000的栅格。

若想提高插值精度可以提高栅格的高度和宽度，当然计算量会变大。这里需要根据具体数据特征和研究需求综合考虑，需要保证插值结果的精度和计算效率。

// 获取预警区域的范围
Geometry cropGeometry = getCropGeometry();
// 根据预警区域获取栅格的范围
Envelope envelope = JTS.toEnvelope(cropGeometry);
double xMin = envelope.getMinX();
double xMax = envelope.getMaxX();
double yMin = envelope.getMinY();
double yMax = envelope.getMaxY();
// 设置栅格 宽度,高度
int height = 1000;
int width = 1000;

模型训练

读取气象站点数据，构建参数，使用克里金插值方法，完成模型训练。

// 构建克里金插值分析
Kriging kriging = new Kriging(Kriging.SPHERICAL_MODEL, 0, 100);
// 获取需要插值的属性值数组
double[] targetValues = new double[data.length];
double[] xList = new double[data.length];
double[] yList = new double[data.length];
for(int i = 0; i < data.length; i++) {
    double[] point = data[i];
    targetValues[i] = point[2];
    xList[i] = point[0];
    yList[i] = point[1];
}
// 变异函数拟合
kriging.train(targetValues, xList, yList);

插值预测、给栅格赋值

遍历栅格，使用克里金插值预测每个栅格的温度值并且赋值。

 // 将数据写入 WritableRaster 对象中
for (int x = 0; x < width; x++) {
    for (int y = 0; y < height; y++) {
        double zbx = xMin + (xMax - xMin) * x / width;
        //raster 是从左上角开始的，所以这里从上往下
        double zby = yMax - (yMax - yMin) * y / height;
        double value  = kriging.predict(zbx,zby);
        raster.setSample(x,  y, 0, value);
    }
}

Geotiff文件保存

将栅格文件保存为geotiff类型文件。

GeoTIFF文件是一种常用的地理信息文件格式，它是标准的TIFF（Tagged Image File Format）格式的扩展，除了存储图像数据以外，还可以包含地理坐标信息和地图投影信息等元数据。因此，GeoTIFF文件可以将图像数据和地理信息完整地存储在一个文件中，方便地图制图、地理信息系统等领域的应用。

// 构建gridCoverage2d
GridCoverageFactory factory = new GridCoverageFactory();
GridCoverage2D coverage = factory.create("contour_polygon", raster,
        new Envelope2D(crs, xMin, yMin, xMax - xMin, yMax - yMin));
// 保存栅格文件
GeoTiffWriter writer = new GeoTiffWriter(outputFile);
writer.write(coverage, null);
writer.dispose();

栅格使用单波段灰度渲染，单波段伪彩色蓝红渲染。

栅格裁剪

如果需要在可视化阶段需要用到栅格文件，可以对栅格文件按掩膜（预测区域）进行裁剪。

public GridCoverage2D crop(GridCoverage2D coverage,  Geometry  geometry) throws IOException {
    CropCoverage cropCoverage = new CropCoverage();
    // 裁剪,返回裁剪后的栅格
    // progressListener 为进度监听器，可以为null
    GridCoverage2D crop = cropCoverage.execute(coverage,geometry,null);
    return crop;
}

栅格发布

通过GeoServer Rest API将处理好的GeoTIFF文件自动化发布到GeoServer地图服务器上，供可视化使用。

public static boolean releaseTiff(String storeName, String fileUrl) throws MalformedURLException, FileNotFoundException {
    URL u = new URL(url);
    //创建一个geoServer rest连接对象
    GeoServerRESTManager manager = new GeoServerRESTManager(u, geoUsername, geoPassword);
    //判断数据存储桶是否存在
    RESTDataStore restStore = manager.getReader().getDatastore(imageWorkspace, storeName);
    //如果不存在就创建一个数据存储，并发布
    if (restStore == null) {
        GSGeoTIFFDatastoreEncoder gsGeoTIFFDatastoreEncoder = new GSGeoTIFFDatastoreEncoder(storeName);
        gsGeoTIFFDatastoreEncoder.setWorkspaceName(imageWorkspace); 
        boolean createStore = manager.getStoreManager().create(imageWorkspace, gsGeoTIFFDatastoreEncoder);
        log.info("create store (TIFF文件创建状态) : " + createStore);
        boolean publish = false;
        publish = manager.getPublisher().publishGeoTIFF(imageWorkspace, storeName, storeName, new File(fileUrl));
        log.info("publish (TIFF文件发布状态) : " + publish);
        if (publish) {
            return true;
        }
    } else {
        log.info("数据存储已经存在了,store:" + storeName);
    }
    return false;
}

制作等值面

确定等值面区间

参考中国气象局网站，设置最低温-40度，最高温40度。每4度为一个区间，构建20个区间。

public static List<Range> tempRangeList = new ArrayList<>();
static {
    // 最高温 40 最低温 -40,每4度一个区间, 参考中国气象局
    for(int i = -40 ; i < 40 ; i = i + 4){
        tempRangeList.add(Range.create(i,true, i+4,false));
    }
}

栅格转面

读取GeoTiff文件，根据等值面区间栅格转成等值面，完成温度等值面专题图。

public SimpleFeatureCollection polygonExtract(GridCoverage2D coverage,List<Range> ranges)  throws Exception{
    PolygonExtractionProcess polygonExtractionProcess = new PolygonExtractionProcess();
    SimpleFeatureCollection polygonFeatures = polygonExtractionProcess.execute(coverage,0,true,null,null,ranges,null);
    return polygonFeatures;
}

数据入库

将温度等值面保存至postgis中。

shapeFileManager.saveToPostGIS(intersectionCollection, 
                               "tmin_polygon_intersection",reportDate);

图层发布

可以选用以下面两个方法任意一种进行图层发布。

（1）通过GeoServer Rest API将处理好的等值面shape文件自动化发布到GeoServer地图服务器上，供可视化使用。

（2）可以将PostGIS中的空间数据表作为源发布到GeoServer中，可视化时候根据reportDate进行筛选，选择想要展示的等值面。

相交分析、生成预警路段

相交分析

将高速公路图层与等值面图层进行相交处理。

可以选用以下面两个方法任意一种进行图层相交处理。

（1）使用PostGIS提供的空间函数，进行图层相交分析。

使用空间数据库的提供的函数，会简化代码。但是需要注意的是，ST_Intersection 函数在 PostGIS 3.0 及以上版本中可以支持多线程并行计算。 如果是大规模数据处理时，并非所有的情况都适合使用多线程计算，如果数据量过大，需要分布式计算的时候，查看数据库集群是否支持。

SELECT
	st_astext(ST_Intersection ( tmin_polygon_intersection.the_geom, highway_polyline.geom )) AS INTERSECTION,
	highway_polyline."NAME",
	tmin_polygon_intersection."value" 
FROM
	tmin_polygon_intersection,
	highway_polyline 
WHERE
	ST_Intersects ( tmin_polygon_intersection.the_geom, highway_polyline.geom ) 
AND 
	tmin_polygon_intersection."reportDate" = '2021-11-21'

（2）采用geotools进行图层相交分析。

此方法会增加代码复杂性。优点是在处理大规模数据的时候，可以配合大数据框架。

// 裁剪
SimpleFeatureCollection intersectionCollection =
    shapeFileManager.intersection(polygonFeatures, china_region);
// 保存至PostGIS
shapeFileManager.
    saveToPostGIS(intersectionCollection, "D://test//intersection.shp");

生成预警路段

可视化

WebGIS

专题图导出

专题图导出的难点在配图，需要花大量时间去调整地图样式，图例等。这里只是展示了简单的配图，对美观度有需求的，需要自己去配置样式。

// 按照value值 做地图渲染
StyleFactory styleFactory = CommonFactoryFinder.getStyleFactory();
FilterFactory2 filterFactory = CommonFactoryFinder.getFilterFactory2();
// 创建边框样式（Stroke）
Stroke stroke = styleFactory.createStroke(
    filterFactory.literal(Color.BLACK), // 边框颜色
    filterFactory.literal(1), // 边框宽度
    filterFactory.literal(1)// 边框透明度
);
FeatureTypeStyle fts = styleFactory.createFeatureTypeStyle();
// 添加rgb值
int[][] rgbs = new int[10][3];
rgbs[0] = new int[]{101,133,143};
rgbs[1] = new int[]{59, 227, 233}; 
// 这里偷懒，只做了2级至11级的样式，实际上应该根据实际情况做更多的样式
for( int i = 2 ; i <= 11 ; i++){
    int[] rgb = rgbs[i-2];
    // 创建填充样式（Fill）
    Fill fill = styleFactory.createFill(
        filterFactory.literal( new Color(rgb[0],rgb[1],rgb[2]))
    );
    Rule rule = styleFactory.createRule();
    Function parseDoubleFunction = filterFactory.function("parseDouble", filterFactory.property("value"));
    rule.setFilter(filterFactory.equals(parseDoubleFunction, filterFactory.literal(i)));
    Symbolizer symbolizer = styleFactory.createPolygonSymbolizer(stroke, fill, "the_geom");
    rule.symbolizers().add(symbolizer);
    fts.rules().add(rule);
}
Style style = styleFactory.createStyle();
style.featureTypeStyles().add(fts);

温度专题图

叠加温度等值线、创建样式文件，确定图片大小，地理信息，导出专题图。

高速公路结冰预警专题图(补充)

预警报告(补充）

将相交分析结果在 GIS 软件中进行可视化，清晰地显示出潜在结冰的高速公路区域。为这些区域发布结冰预警，提醒驾驶员注意行车安全。