利用kriging.js和turf.js生成等值线图的geojson

前言

最近项目有一个需求,要求使用前端生成等值线图数据,并且因为甲方使用的自己的三维地图显示sdk,所以,我不得不离开舒适的cesium去研究纯js的geojson生成,总而言之对于我这种地理信息处理半桶水还是有点难度的,但是好在最后还是没出啥问题的搞出来了,于是记录一下

获取geojson

总体上来讲,分为三步

• 第一步,使用kriging.train生成variogram对象
• 第二步,使用kriging.grid生成栅格化数据,推导每个栅格的值是多少,并将值在一个区间内的栅格放到一个多边形内

1. 拿到的数据

我们手头能拿到的数据有点位数据,即经纬度加上数值,大概是这样

[    {        lat: 'xxx',        lng: 'xxx',        value: 'xxx'    },    ...]

还有地图的边界数据,即渲染等值线的总区域的边界线数据,用于生成合适大小的矩形等值线图和最后裁剪

2. 获取kriging.train()的参数

有了这两个数据之后,我们首先处理点位数据,krigingjs的train方法要求传入经度数组

这个是前端基本功,通过for循环分别将value,lat,lng的值 push到 values/lats/lngs 三个数组中,代码我就不写了

3. 生成variogram对象

没啥好说的,直接照着写

    // 参数照着写就行,我也不懂
    let params = {
        //model可选'gaussian(高斯模型)','spherical(球形模型)',exponential(指数模型) 
        krigingModel: 'exponential',
        krigingSigma2: 0,
        krigingAlpha: 100,
        canvasAlpha: 0.75,//canvas图层透明度
    };
    let variogram = kriging.train(
        values,
        lons,
        lats,
        params.krigingModel,
        params.krigingSigma2,
        params.krigingAlpha
    );
    console.log('variogram', variogram);

4. 获取kriging.grid()的参数,并生成grid对象

需要polygons/variogram/width三个参数

• polygons代表生成的等值线图的矩形边界
• variogram就是处理过的点位数据
• width是栅格的宽度

width看着填就行,也有算出来的,但我觉得没啥必要,所以现在只需要获取polygons,使用turf.bbox()方法获取

    import boundaries from '@/assets/js/haishu1.json'
    
    let polygons = [];
    const extent = turf.bbox(boundaries)
    polygons.push([
        [extent[0], extent[1]], [extent[0], extent[3]],
        [extent[2], extent[3]], [extent[2], extent[1]]
    ]);
    
    let grid = kriging.fastGrid(polygons, variogram, 0.0015);

5. 利用网格计算点集

计算栅格的点对应的经纬度和值

const gridFeatureCollection = function (grid, xlim, ylim) {
    var range = grid.zlim[1] - grid.zlim[0];
    var i, j, x, y, z;
    var point
    var n = grid.length;//列数
    var m = grid[0].length;//行数
    var pointArray = [];
    for (i = 0; i < n ; i++)
        for (j = 0; j < m ; j++) {
            x = (i) * grid.width + grid.xlim[0];
            y = (j) * grid.width + grid.ylim[0];
            z = +grid[i][j];
            const isNaN = Number.isNaN(z)
            if(isNaN) {
                // 有时候z可能是NaN,应该是数据本身有问题,这里我看缺的部分周围的值,给了一个中间值
                z = 1.2494884579461856
            }
            
            point = turf.point([x, y], { value: z })
            pointArray.push(point);
        }
    return pointArray;
}

let fc = gridFeatureCollection(grid);

6.利用turf.isobands()生成等值面的多边形

需要的参数 :

• collection : 点集,生成等值线图的数据支撑
• scaleList : 等值面颜色
• 其他设置,下面使用的有
  • zProperty: 取features的properties中的哪个属性为值
  • commonProperties: 通用设置,fill-opacity为透明度设置
  • breaksProperties: 设置等值的区间,同一区间内的值视为等值

a. collection的获取

利用turf画等值面,collection参数要处理成geojson类型,将点集处理一下

    var collection = turf.featureCollection(fc);
    console.log('collection', collection);

b. scaleList的获取

项目要求有常规展示和随值展示两种,常规展示为固定区间,随值展示会根据数据中值的大小获取区间

总而言之,固定值叫产品给你,随值问产品需要多少个刻度然后最大之间最小值除一下push到数组里就行

// 结果大概是这种
let scaleList = [-10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45]

c. 其他设置

    var isobands_options = {
        zProperty: "value",
        commonProperties: {
            "fill-opacity": 0.8  //设置通用属性
        },
        breaksProperties: [
                    {fill:'#40a4c0'},
                    {fill:'#6ec2bf'},
                    {fill:'#92c1a5'},
                    {fill:'#bab37c'},
                    {fill:'#c3af71'},
                    {fill:'#f29d40'},
                    {fill:'#ef6c31'},
                    {fill:'#ee6932'},
                    {fill:'#e65231'},
                    {fill:'#af444b'},
                ]  // 设置每个区间的颜色
    };

d. 生成等值线图的geojson,并优化

    var isobands = turf.isobands(
        collection,
        opt.scaleList,  // 设置值的颜色区间
        isobands_options
    );
    console.log('isobands', isobands);
    //按照面积对图层进行排序，规避turf的一个bug
    function sortArea(a,b) {
        return turf.area(b)-turf.area(a);
    }
    isobands.features.sort(sortArea)
    // turf.flatten 扁平化数据,减少嵌套层级
    boundaries = turf.flatten(boundaries);  
    isobands = turf.flatten(isobands);

e. 对于超出边界的区域进行裁剪

var features = [];
    
    isobands.features.forEach(function (layer1) {
      boundaries.features.forEach(function (layer2) {
        let intersection = null;
        try {
          intersection = turf.intersect(layer1, layer2);
        }
         catch (e) {
          layer1 = turf.buffer(layer1, 0);
          intersection = turf.intersect(layer1, layer2);
        }
        if (intersection != null) {
          intersection.properties = layer1.properties;
          intersection.id = Math.random() * 100000;
          features.push(intersection);
        }
      });
    });
    
    var intersection = turf.featureCollection(features);

得到的 intersection 就是我们要的最终的等值线的geojson了,渲染的话,直接通过cesium的api加载就行,如果你和我一样用的sdk没有这么方便的功能,那么你需要根据 intersection 中的多边形数组一个一个渲染,还要注意有孔的和多面的多边形,总之对于不熟悉geojson中属性带变得含义的人来说还是很麻烦的,以防万一,我也把我处理的函数写下来

顺便放一下geojson属性的解释# GeoJSON 对象

注 : earth参数是地图sdk的实例,不需要太关注,data是上面做好的等值线图的geojson

export function main(earth, data) {
    for(let i = 0; i<data.features.length; i++) {

        const item = data.features[i];
        item.geometry.coordinates.forEach(item1 => {
            let {points,interPoints} = getPoint(item1)
            console.log('points =======> ',points);
            const opt = {...item.properties}
            createElementPolygon(earth, {points,interPoints}, opt)
        })
    }
}

function getPoint(item1){
    let points = [],interPoints=[]
    //    面
    let point = item1[0];
    // 有孔
    if(point.length > 2){
        item1.forEach((item2,index) => {
            if(index === 0){
                item2.forEach(item3 => {
                    let [x,y] = item3
                    points.push(StampGis.Cartesian3.fromDegrees(x, y, 0))
                })
            }else{
                let interPoint = []
                item2.forEach(item3 => {
                    let [x,y] = item3
                    interPoint.push(StampGis.Cartesian3.fromDegrees(x, y, 0))
                })
                interPoints.push(interPoint)
            }
        })

    }
    else{
        // 无孔
        item1.forEach(item2 => {
            let [x,y] = item2
            points.push(StampGis.Cartesian3.fromDegrees(x, y, 0))
        })
    }
    return {points,interPoints}
}

BUG

在使用的过程中,我发现数据多数处于一个区间内的时候,该区间可能会生成大量空隙,能看到turf.isobands()拿到的geojson中feature.geometry下coordinates是空数组,暂时不知道是什么原因,所以也没解决,所以有大佬看到这篇文章有兴趣能指点一下真的拜谢了!

目前我使用的方法是通过计算所有值的平均值所在的区间,将该区间的颜色作为底图添加在等值线图的最底层,算是没办法的办法,方法我写在下面了

function drawIso() {
   //...之前的代码需要将colorList单独拿出来定义一个变量
       const scaleList = this.getScaleList(max,min)
            const colorList = [
                    {fill:'#04007a'},
                    {fill:'#121aac'},
                    {fill:'#253fd1'},
                    {fill:'#2569dc'},
                    {fill:'#3084ed'},
                    {fill:'#58b9f4'},
                    {fill:'#91d8ff'},
                    {fill:'#b8eff1'},
                ]
            console.log('scaleList', scaleList);
            const isoData = isoLine({
                colorList,
                scaleList,
                krikingOpt: {
                    lons,lats,values
                },
            })
            console.log('isoData',isoData);
            let colorIndex
            const average = values.reduce((pre, item) => {
                return pre + +item
            },0) / values.length
            scaleList.forEach((item,index) => {
                if(+item >= average && +(scaleList[index+1]) <= max) colorIndex = index
            })
            main(this.stampAPI.usearth, isoData, colorList[colorIndex].fill)
}

在之前写的main函数中最开始调用isoBaseLayer()函数

export function main(earth, data) {
    isoBaseLayer(earth, isoBaseLayerFill)
    // ....剩下的
}

最后是添加底图的方法

function isoBaseLayer(earth,fill) {
    if(!fill) return
    const data = boundaries1.features[0].geometry.coordinates[0].map(item => {
        return StampGis.Cartesian3.fromDegrees(item[0], item[1], 50)
    })
    // createElementPolygon就是接收多边形的笛卡尔3数组,创建多边形,使用你自己的创建多边形函数
    createElementPolygon(earth, {points:data}, {
        'fill-opacity': 0.5, fill: fill
    })
}