OpenLayers结合Turf实现空间运算

1. 引言

空间运算利用几何函数来接收输入的空间数据，对其进行分析，然后生成输出数据，输出数据为针对输入数据执行分析的派生结果。

可从空间运算中获得的派生数据包括：

作为输入要素周围缓冲区的面
作为对几何集合执行分析的结果的单个要素
作为比较结果以确定不与其他要素位于同一物理空间的要素部分的单个要素
作为比较结果以查找与其他要素的物理空间相交的要素部分的单个要素
由彼此不位于同一物理空间的输入要素部分组成的多部分 (multipart) 要素
作为两个几何的并集的要素

参考文档：空间运算—ArcMap | 文档 (arcgis.com)

Turf.js是MapBox公司研发的基于浏览器端的空间分析库，它使用JavaScript进行编写，通过npm进行包管理。值得一提的是，良好的模块化设计使其不仅能够作用于浏览器端、还可通过Node.js在服务端使用。Turf 原生支持 GeoJSON 矢量数据。GeoJSON 的优点是结构简单，并且得到了所有网页地图API的支持；但 GeoJSON 不支持空间索引，这个缺点可能会限制 Turf 处理大型文件的能力效率。其适用于轻量级（数据轻量而非功能轻量）的WebGIS应用

参考文献：Turf.js—让你在浏览器上实现地理分析 - 掘金 (juejin.cn)

turf.js官网：Turf.js | Advanced Geospatial Analysis (turfjs.org)

中文站点：Turf.js中文网 (fenxianglu.cn)

Turf.js Github地址：Turfjs/turf: A modular geospatial engine written in JavaScript (github.com)

空间运算可谓是空间分析的基础，Turf.js提供了大量的空间分析功能，包含了空间运算功能，本文参考OpenLayers与Turf集成的官方示例，使用示例数据和原生JavaScript，进行求交运算与缓冲区运算，并进行可视化

OpenLayers与Turf集成示例：turf.js (openlayers.org)

2.初始化地图

引入OpenLayers和Turf.js的CDN：

<link rel="stylesheet"  href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/css/ol.css"  type="text/css">
<script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/build/ol.js"></script>
<script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script>

注意：

OpenLayers与Turf集成示例使用的Turf版本是2.0.0，版本偏老，与目前版本（6.5.0）部分API不兼容，不建议使用

构建网页基本内容，添加地图容器：

<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"><meta http-equiv="X-UA-Compatible" content="IE=edge"><meta name="viewport" content="width=device-width, initial-scale=1.0"><title>Document</title><link rel="stylesheet"href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/css/ol.css" type="text/css"><script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/build/ol.js"></script><script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script><style>html,body,#map {height: 100%;}</style>
</head><body><div id="map"></div></body></html>

添加地图并加载示例矢量数据:

<script>fetch('https://openlayers.org/en/latest/examples/data/geojson/roads-seoul.geojson').then(response => response.json()).then(function (json){var vectorSource = new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(json, {featureProjection: 'EPSG:3857'})});var vectorLayer = new ol.layer.Vector({source: vectorSource});var map = new ol.Map({layers: [new ol.layer.Tile({source: new ol.source.OSM()}),vectorLayer],target: 'map',view: new ol.View({center: ol.proj.fromLonLat([126.980366, 37.52654]),zoom: 15})});})
</script>

完整代码如下:

<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"><meta http-equiv="X-UA-Compatible" content="IE=edge"><meta name="viewport" content="width=device-width, initial-scale=1.0"><title>Document</title><link rel="stylesheet"href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/css/ol.css" type="text/css"><script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/build/ol.js"></script><script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script><style>html,body,#map {height: 100%;}</style>
</head><body><div id="map"></div><script>fetch('https://openlayers.org/en/latest/examples/data/geojson/roads-seoul.geojson').then(response => response.json()).then(function (json){var vectorSource = new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(json, {featureProjection: 'EPSG:3857'})});var vectorLayer = new ol.layer.Vector({source: vectorSource});var map = new ol.Map({layers: [new ol.layer.Tile({source: new ol.source.OSM()}),vectorLayer],target: 'map',view: new ol.View({center: ol.proj.fromLonLat([126.980366, 37.52654]),zoom: 15})});})</script></body></html>

初始化的地图:

3. 缓冲区计算

参考官方文档:Turf.js | Advanced geospatial analysis (turfjs.org)

步骤实质就是:

将Openlayers的Features对象转换为Turf支持的GeoJSON对象
调用Turf.buff()函数得到运算结果
OpenLayers读取运算后的GeoJSON数据加载到图层中

缓存区计算的核心代码如下:

const features = vectorSource.getFeatures();
const turfLines = (new ol.format.GeoJSON()).writeFeaturesObject(features, {featureProjection: 'EPSG:3857'
});
var buffered = turf.buffer(turfLines, 25, { units: 'meters' });
const bufferedLayer = new ol.layer.Vector({source: new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(buffered, {featureProjection: 'EPSG:3857'})})
})
map.addLayer(bufferedLayer);

注意:

从OpenLayers的Feature对象转换为Turf支持的GeoJSON对象时务必指定坐标系

这一步的完成代码如下:

<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"><meta http-equiv="X-UA-Compatible" content="IE=edge"><meta name="viewport" content="width=device-width, initial-scale=1.0"><title>Document</title><link rel="stylesheet"href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/css/ol.css" type="text/css"><script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/build/ol.js"></script><script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script><!-- <script src="https://api.tiles.mapbox.com/mapbox.js/plugins/turf/v2.0.0/turf.min.js"></script> --><style>html,body,#map {height: 100%;}</style>
</head><body><div id="map"></div><script>var map = new ol.Map({layers: [new ol.layer.Tile({source: new ol.source.OSM()})],target: 'map',view: new ol.View({center: ol.proj.fromLonLat([126.980366, 37.52654]),zoom: 15})});fetch('https://openlayers.org/en/latest/examples/data/geojson/roads-seoul.geojson').then(response => response.json()).then(function (json) {var vectorSource = new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(json, {featureProjection: 'EPSG:3857'})});var vectorLayer = new ol.layer.Vector({source: vectorSource,style: new ol.style.Style({stroke: new ol.style.Stroke({color: '#0000ff',width: 2})})});map.addLayer(vectorLayer);const features = vectorSource.getFeatures();const turfLines = (new ol.format.GeoJSON()).writeFeaturesObject(features, {featureProjection: 'EPSG:3857'});var buffered = turf.buffer(turfLines, 25, { units: 'meters' });const bufferedLayer = new ol.layer.Vector({source: new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(buffered, {featureProjection: 'EPSG:3857'})})})map.addLayer(bufferedLayer);})</script>
</body></html>

结果图如下:

4. 相交计算

参考官方文档:Turf.js | Advanced geospatial analysis (turfjs.org)

步骤与缓冲区计算类似:

将Openlayers的Features对象转换为Turf支持的GeoJSON对象
调用Turf.lineIntersect()函数得到运算结果
OpenLayers读取运算后的GeoJSON数据加载到图层中

本文所使用的数据是LineString,所以此处的相交计算的是线段相交

相交运算的核心代码如下:

const features = vectorSource.getFeatures();
const turfLines = (new ol.format.GeoJSON()).writeFeaturesObject(features, {featureProjection: 'EPSG:3857'
});const intersectionSource = new ol.source.Vector()
for (let i = 0; i < turfLines.features.length; i++) {for (let j = i + 1; j < turfLines.features.length; j++) {const intersections = turf.lineIntersect(turfLines.features[i], turfLines.features[j]);if (intersections) {if (intersections.features.length > 0) {intersectionSource.addFeatures(new ol.format.GeoJSON().readFeatures(intersections, {featureProjection: 'EPSG:3857'}))}}}
}const intersectionLayer = new ol.layer.Vector({source: intersectionSource,style: new ol.style.Style({stroke: new ol.style.Stroke({color: '#ff0000',width: 2}),fill: new ol.style.Fill({color: '#ff0000'}),image: new ol.style.Circle({radius: 5,fill: new ol.style.Fill({color: '#ff0000'})})})
})
map.addLayer(intersectionLayer);

注意:

计算结果为点集合(Point),最好设置一定的显示样式,不然难以看见

这一步的完成代码如下:

<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"><meta http-equiv="X-UA-Compatible" content="IE=edge"><meta name="viewport" content="width=device-width, initial-scale=1.0"><title>Document</title><link rel="stylesheet"href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/css/ol.css" type="text/css"><script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/build/ol.js"></script><script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script><!-- <script src="https://api.tiles.mapbox.com/mapbox.js/plugins/turf/v2.0.0/turf.min.js"></script> --><style>html,body,#map {height: 100%;}</style>
</head><body><div id="map"></div><script>var map = new ol.Map({layers: [new ol.layer.Tile({source: new ol.source.OSM()})],target: 'map',view: new ol.View({center: ol.proj.fromLonLat([126.980366, 37.52654]),zoom: 15})});fetch('https://openlayers.org/en/latest/examples/data/geojson/roads-seoul.geojson').then(response => response.json()).then(function (json) {var vectorSource = new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(json, {featureProjection: 'EPSG:3857'})});var vectorLayer = new ol.layer.Vector({source: vectorSource,style: new ol.style.Style({stroke: new ol.style.Stroke({color: '#0000ff',width: 2})})});map.addLayer(vectorLayer);const features = vectorSource.getFeatures();const turfLines = (new ol.format.GeoJSON()).writeFeaturesObject(features, {featureProjection: 'EPSG:3857'});const intersectionSource = new ol.source.Vector()for (let i = 0; i < turfLines.features.length; i++) {for (let j = i + 1; j < turfLines.features.length; j++) {const intersections = turf.lineIntersect(turfLines.features[i], turfLines.features[j]);if (intersections) {// console.log(intersections)if (intersections.features.length > 0) {intersectionSource.addFeatures(new ol.format.GeoJSON().readFeatures(intersections, {featureProjection: 'EPSG:3857'}))// console.log(intersectionSource.getFeatures())}}}}// console.log(intersectionSource)const intersectionLayer = new ol.layer.Vector({source: intersectionSource,style: new ol.style.Style({stroke: new ol.style.Stroke({color: '#ff0000',width: 2}),fill: new ol.style.Fill({color: '#ff0000'}),image: new ol.style.Circle({radius: 5,fill: new ol.style.Fill({color: '#ff0000'})})})})map.addLayer(intersectionLayer);})</script>
</body></html>

结果图如下:

5. 整合

将两个运算结果结合在一起,完整代码如下:

<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"><meta http-equiv="X-UA-Compatible" content="IE=edge"><meta name="viewport" content="width=device-width, initial-scale=1.0"><title>Document</title><link rel="stylesheet"href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/css/ol.css" type="text/css"><script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.14.1/build/ol.js"></script><script src='https://unpkg.com/@turf/turf@6/turf.min.js'></script><!-- <script src="https://api.tiles.mapbox.com/mapbox.js/plugins/turf/v2.0.0/turf.min.js"></script> --><style>html,body,#map {height: 100%;}</style>
</head><body><div id="map"></div><script>var map = new ol.Map({layers: [new ol.layer.Tile({source: new ol.source.OSM()})],target: 'map',view: new ol.View({center: ol.proj.fromLonLat([126.980366, 37.52654]),zoom: 15})});fetch('https://openlayers.org/en/latest/examples/data/geojson/roads-seoul.geojson').then(response => response.json()).then(function (json) {var vectorSource = new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(json, {featureProjection: 'EPSG:3857'})});var vectorLayer = new ol.layer.Vector({source: vectorSource,style: new ol.style.Style({stroke: new ol.style.Stroke({color: '#0000ff',width: 2})})});map.addLayer(vectorLayer);const features = vectorSource.getFeatures();const turfLines = (new ol.format.GeoJSON()).writeFeaturesObject(features, {featureProjection: 'EPSG:3857'});// console.log(turfLines)// const buffered = turf.buffer(turfLines, 25);var buffered = turf.buffer(turfLines, 25, { units: 'meters' });const bufferedLayer = new ol.layer.Vector({source: new ol.source.Vector({features: (new ol.format.GeoJSON()).readFeatures(buffered, {featureProjection: 'EPSG:3857'})})})map.addLayer(bufferedLayer);const intersectionSource = new ol.source.Vector()for (let i = 0; i < turfLines.features.length; i++) {for (let j = i + 1; j < turfLines.features.length; j++) {const intersections = turf.lineIntersect(turfLines.features[i], turfLines.features[j]);if (intersections) {// console.log(intersections)if (intersections.features.length > 0) {intersectionSource.addFeatures(new ol.format.GeoJSON().readFeatures(intersections, {featureProjection: 'EPSG:3857'}))// console.log(intersectionSource.getFeatures())}}}}// console.log(intersectionSource)const intersectionLayer = new ol.layer.Vector({source: intersectionSource,style: new ol.style.Style({stroke: new ol.style.Stroke({color: '#ff0000',width: 2}),fill: new ol.style.Fill({color: '#ff0000'}),image: new ol.style.Circle({radius: 5,fill: new ol.style.Fill({color: '#ff0000'})})})})map.addLayer(intersectionLayer);})</script>
</body></html>

结果图如下:

6. 参考资料

[1]空间运算—ArcMap | 文档 (arcgis.com)

[2]Turf.js—让你在浏览器上实现地理分析 - 掘金 (juejin.cn)

[3]Turf.js | Advanced Geospatial Analysis (turfjs.org)

[4]Turf.js中文网 (fenxianglu.cn)

[5]Turfjs/turf: A modular geospatial engine written in JavaScript (github.com)

[6]使用 Fetch - Web API 接口参考 | MDN (mozilla.org)

[7]turf.js (openlayers.org)

[8]openlayers+turf.js实现缓冲区的绘制_gis_SSS的博客-CSDN博客_openlayers缓冲区分析

[9]OpenLayers v6.14.1 API - Class: EsriJSON