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

作者:yangjunlin

  在上一篇文章中我们已经使用了像素法实现克里金插值的方式，但是问题也就随之抛了出来。1.第一点，在反距离权重插值的时候，因为处理的数据量大会直接导致主线程卡，导致用户体验不好，2.第二点，分析过程中 zoom级别改变，坐标转像素的值的计算方式也会发生变化，这样就会导致(坐标转像素坐标，像素坐标转坐标)会出现错误结果。
  我们先分析一下第一个问题，主线程卡顿的问题其实很好解决， Web Worker 可以创建一个独立于主线程运行的子线程。可以将一些【可能会阻塞主线程的操作】，丢在 Worker 里去单独执行。因此我们可以直接将卡住进程的反距离权重插值分析方法直接抛到webworker里面即可
  再看下第二个问题，既然像素没法用，那么我们就自己搞‘像素’ 也就是自己画格子。其实坐标到像素仔细来看可以认为是一种投影，那么我们是不是也能自己搞一套投影了。从当前我的数据的坐标系来看，我用的mapboxgl3857的web墨卡托投影的坐标，全球的比例级别可直接认为是1:1,中国地区是没有被压扁的。我的思路是这样的我们也可以画一个1:1的范围是正方形的图形，然后对里面的区域再进行[m,m]的小网格划分。3857web墨卡托与该正方形的范围基本上就是成一个k比例的，也就是说，我们不用复杂的投影计算，只需要简单的算出每个单位的经纬度所占的小格子数就可以，这样就大大降低了难度。

【第二问题的解决】下面我将带着大家先搞第二个问题，用自己创建网格的方式来做前端克里金插值分析

• 1. 首先我们先搞正方形的范围，我采取的方案是将分析范围矩形变为正方形，这个很简单直接用矩形的长边作为正方形的边长即可，然后拿到改边的长度

 let width=bbox[2]-bbox[0];
    let height=bbox[3]-bbox[1];
    length=width-height>0?width:height
    return length;
1
2
3
4

• 2. 确定好网格的大小，我这里使用是的1000*1000

//形成1000*1000的网格。
    let cw=1000,ch=1000;

1
2
3

• 3. 拿到网格和目前坐标系单位的比例关系

let k=cw/lenght; 
1

• 4. 将气象点坐标，放在几何范围左下为原点的位置上，然后进行所谓的‘投影’–将坐标系单位，转为网格位置，这样我们就能将我们画好的网格和网格坐标系下的分析点的位置关联起来，也就是说后面可以进行反距离权重分析了，将网格上的每个小格子上都填上值

 //将几何范围放到左下为原点的位置上
       let x=data[i].x-xoffset;
       let y=data[i].y-yoffset;
 //进行投影，将局部 
       let posx=x*k;
       let posy=y*k;
1
2
3
4
5
6

• 5.进行反距离权重，也就是填满小格子的操作(前文有描述，这里不在写)
• 6. 然后利用d3-couter得到等值线，特别需要注意的是这里得到的只是格网单位上的值，而不是真实地理位置上的值，我们需要进行一个’反投影’并再加载几何范围的最小x和最小y即可

let x= coordinates[i][j][m][0]*1/k+box[0]
let y=coordinates[i][j][m][1]*1/k+box[1]
coordinates[i][j][m] = [x, y]
1
2
3

• 7. 剩下的就只是相交和渲染了与上文所描述的一致我们现在看下效果，结果发现有偏移，这是为啥呢
     
     其实很简单，就是格子画的不够小，导致分析结果有误差，现在只要我们在细化下网格搞一个5000*5000的，得到的结果非常的完美
     
     
     与像素化的分析结果基本是一致的。但是现在这个我们已经可以控制你划分的格网大小了

【第一个问题解决】这个问题就比较简单了，直接搞个子线程将反距离权重交给他来做即可直接上代码

• 1. 首先创建worker.js,并进行监听，监听主线程传过来进行反距离分析的数据data，然后再将结果数据返回给主线程

//Worker.js
self.onmessage = function ({data}) {
 const resdata=IDW(data['data'],data['cw'],data['ch'])
 self.postMessage({ resdata })
}
self.onerror = function (err) {
 throw new Error(err)
}
1
2
3
4
5
6
7
8

//主进程：
//webworker
const worker = new Worker('./worker.js', {
   name: 'my-worker'
})
//主线程向子线程发数据
worker.postMessage({
   data: pxdata.data,
   cw: cw,
   ch: ch,
   contours
})
// 监听到子线程返回的数据
worker.onmessage = ({
   data
}) => {
   const idwdata = data['resdata'];}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

然后将拿到结果在进行d3分析，解析渲染即可，这样就解决因为反距离权重分析的时间长导致主进程卡在而导致无法正常使用地图的问题。

后言

当然这个只是我的一个demo，代码肯定还有很多不完善的地方，这两篇文章只是提供给大家做前端插值分析的一个小方向，如果大家有更多的更优的解决方案，欢迎大家在评论里讨论