前端实现克里金插值分析(一)

作者:yangjunlin

最近不少小伙伴问我怎么搞前端插值分析，我在github上查找了一些资料，目前最常用的方式是webgis框架+idw（反距离权重算法）+d3-contour的方式实现，这种方式是比较简单同时基本能满足一般的气象分析，比如局部的平均降水发布图等，下来我将详细的讲述具体的实现方式

首先我们需要知道什么是d3-contour,具体的使用方式是什么，可以参考https://www.wenjiangs.com/doc/d3-contour-2这个文档，阅读之后我们发现他是一个类库，

以我的角度理解是将任意一个矩形，将其分成若干个单位相同的正方形格子，并且记录下每个正方形中的value值（非常类似栅格图像以像素为单位进行分割）。我红框圈起来的就非常能体现，一个256*256的矩形盒子以0.5为单位的正方形进行切割，每个正方形有个值，然后调用方法，只需要传入thresholds阈值返回的并返回每个阈值对应的等值线坐标数组。看到这里是不是有点想法了，我们只要想方法将我们分析的区域搞成矩形，并同时将其分割成若干个单位相同的正方形，然后将采样点和这个‘单位’正方形关联起来。既然提到了单位正方形，就非常类似咋们的像元了，因此废话不多说直接上代码

首先准备数据

全国面数据 2.2011年各地区气象站点数据（属性值中包含降雨量）

具体实现步骤

• 1. 刚才咱们说需要把分析区域搞为矩形区域，能最快能想到的方式是利用turf的bbox拿到矩形范围，当然比如openlayer等都自带可以获得几何矩形范围的方法。因为我使用的iclient for mapboxgl 就行开发的，这个webgis库并没带有该方法，所以我采取的turf获取

import bbox from '@turf/bbox'
let box = bbox(region);

• 2. 下一步就是切割了,经过查找资料发现最快捷的方式其实就是将矩形几何范围转为几何像素的范围，该范围内的单位刚好就是一个像素，像素可以理解为就是正方形的，并且后面站点的坐标也能直接转为像素，也方便站点值与切割的小矩形(像元)进行关联。因此我们现在只需要将得到的几何范围(左上右下)转为像素坐标得到一个像素范围，其实每个单位代表一个像元即上文提到的’小矩形’（当然也有弊端，zoom级别不同，导致单位像元的大小不一致的情况）

let p = map.project(pos);

• 3. 然后我们就能得到上面讲的size(n,m)就是像元范围的上下的差值以及左右的插值

let cw = pxregion.xmax - pxregion.xmin;
let ch = pxregion.ymax - pxregion.ymin;

• 4. 现在就要给像素范围里面的每个小正方形(像素)填值，之前我们提到过可以直接将气象点的坐标转为像素坐标，但是有个注意点，如果要将气象点的值与像素范围上的像素关联起来(类似空间上的传递赋值)是需要坐标一致的，现在我们拿到的像素范围的坐标是以(像素范围的最小x，像素范围最小的y)作为左下的原点，那么现在就需要将气象点的像素坐标值，减去分析像素范围的最小x和像素范围最小的y,来得到统一坐标系下的站点的像素位置和值

let p = map.project([data[i].x, data[i].y])
_data.push({
                   x: p['x'] - pxregion.xmin,
                   y: p['y'] - pxregion.ymin,
                   value: _value
               })

• 5.现在我们已经拿到了n,m分析范围，和n,m上某些像元的值，下一步需要做的就是利用反距离权重来得到，像素范围内的每个像素的值，具体的实现方式网上有很多算法，我取了一个常用的

        const dlen = d.length;
       let matrixData = new Array(width * height);
       let idwcount = 0
       for (let i = 0, k1 = 0; i < height; i++) {
           for (let j = 0; j < width; j++, k1++) {
               let sum0 = 0,
                   sum1 = 0;
               for (let k = dlen - 1; k >= 0; k--) {
                   const dk = d[k];
                   const dis = Math.pow((i - dk.y), 2) + Math.pow((j - dk.x), 2);
                   sum0 += dk.value / dis;
                   sum1 += 1 / dis;
                   idwcount++
               }
               matrixData[k1] = sum0 / sum1
           }
       }
return matrixData;

• 5. 反距离插值以后我们就能得到像素范围内的每个像元的值，然后利用d3-contour得到等值区域。[20, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 11000, 12000, 13000, 14000, 15000, 16000, 17000, 18000, 19000, 20000]，这个是分析的阈值，可根据你们直接想分段的情况进行填写

contours.thresholds([20, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 11000, 12000, 13000, 14000, 15000, 16000, 17000, 18000, 19000, 20000])(idwdata).forEach((fakegeometry) => {})

• 6. 现在返回的是不同阈值下的假的几何坐标和对应的阈值，为什么说假的呢，因为现在的这个坐标还是像素坐标，他不是正真的几何坐标，因此还需要转换将像素坐标转为几何坐标，记得要将以前为了统一坐标而减去的分析像素范围的最小x和像素范围最小的y加上，加上以后再进行转换

let data = map.unproject([coordinates[i][j][m][0] + pxregion.xmin, coordinates[i][j][m][1] + pxregion.ymin])
coordinates[i][j][m] = [data.lng, data.lat]

• 7. 得到最终geometry和对应的阈值以后很快的能将其构造成feature，但是还没结束，现在得到的分析结果依然是个矩形范围，我们需要得到拿来分析的几何的范围，因此还需要将分析出的结果与分析几何进行相交操作，我这边依然用的turf来进行相交

import intersect from '@turf/intersect' 
dataset.features.map(item => {
       let resfea = intersect(item, region.features[0])
       resfea['properties']['value'] = item['properties']['value']
       datasetres.features.push(resfea)
   })

• 8. 最后利用mapboxgl进行要素渲染，我们来看下结果，并且分析出来的结果与采样点对比基本一致
     

后言

最后，这个才做好了一半，为什么说只做好了一半呢，因为有两个大原因1.分析的速度很慢，将主线程卡住了，导致整个界面都没法动，用户体验不好只是不要卡到主进程，2.分析的快慢和你所展示地图的zoom有关系，并且如果当你放的越大的时候，还没法分析出结果或者结果错误(这是因为我们使用的webgis框架自带的坐标转像素坐标导致，这个的结果和当前zoom有直接关系)。我们后面应该采取何种方式避免zoom级别导致的问题并做性能上的优化呢？大家可以关注我的下篇文章，前端克里金插值分析的优化和提升。下篇文章中我将放出全部完整代码及数据(像素渲染和网格渲染都将分别放出)