javascript反距离idw空间插值

标签： Gis 分类： Gis / Javascript 创建时间：2019-11-09 06:12:28 更新时间：2024-11-15 10:49:44

IDW(Inverse Distance Weighted)反距离权重插值，也可以称为距离倒数乘方法。是指距离倒数乘方格网化方法是一个加权平均插值法，可以进行确切的或者圆滑的方式插值。

参考网友 Liquid泉 提供的方法函数和思路，加上我自己的改造，最终实现了前端插值与运算。效果图：

基本的思路是这样的，我获取了点图层，然后根据点图层的属性进行计算插值，然后添加了一个source为ImageCanvas的Image图层，在图层中返回了一个canvas，使用的是d3以及d3-contour进行了等值线的绘制。直接上代码，有点多。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
<!DOCTYPE html>
<html>
<head>
    <title>色斑图测试</title>
    <meta charset="utf-8"/>
    <link rel="stylesheet" href="../lib/ol/ol.css" type="text/css">
    <script src="../lib/ol/ol.js"></script>
    <script  src="https://code.jquery.com/jquery-3.4.1.min.js"  integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script>
    <script src="lib/d3.v4.min.js"></script>
    <script src="lib/d3-contour.v1.min.js"></script>
    <style>
        body,html{
            width:100%;
            margin:0;
            padding:0;
            height: 100%;
            position:absolute;
        }
        .map{
            width:100%;
            height:100%;
            position:absolute;
        }
        .save{
            position: absolute;
            left: 200px;
            z-index: 10;
        }
    </style>
</head>
<body>
    <div class="map" id="map"></div>
    <script src="../lib/lodash.min.js"></script>
    <script src="../js/TMap.js"></script>
    <script>
        var mainMap=new TMap('map',[120.349,30.34],16);
        var vectorSource=null;
        var vector=null;
        var Regioncoord=[];


        //请求数据
        $.ajax({
            url:"js/station.json"
            ,success:function(data){
                var geojsonFormat=new ol.format.GeoJSON();
                var features=geojsonFormat.readFeatures(data);

                //绘制图层
                vectorSource = new ol.source.Vector({
                    name:"station_source",
                    format: new ol.format.GeoJSON(),
                });
                vectorSource.addFeatures(features);
                vector = new ol.layer.Vector({
                    source: vectorSource,
                    name:"station_layer",
                    zIndex:15
                });
                mainMap.addLayer(vector);
                var tempExtent=vectorSource.getExtent();
                Regioncoord.push([tempExtent[0],tempExtent[1]]);
                Regioncoord.push([tempExtent[2],tempExtent[1]]);
                Regioncoord.push([tempExtent[2],tempExtent[3]]);
                Regioncoord.push([tempExtent[0],tempExtent[3]]);

                var count=features.length;
                var pointArray=[];
                for(var i=count-1;i>=0;i--){
                    var feature=features[i];
                    var properties=feature.getProperties();
                    var coord=feature.getGeometry().getCoordinates();

                    var point={};
                    point.Lon=coord[0];
                    point.Lat=coord[1];
                    point.TEM=properties.id;

                    pointArray.push(point);

                }
                // 初始化色斑图
                CanvasLayerInit(pointArray);
            }
            ,error:function(error){
                console.log(error);
            }
        })


        // 色标卡
        var colors = [
                { max: -30, color: '#020c64' },
                { min: -30, max: -28, color: '#071e78' },
                { min: -28, max: -26, color: '#11318b' },
                { min: -26, max: -24, color: '#1b449f' },
                { min: -24, max: -22, color: '#2657b3' },
                { min: -22, max: -20, color: '#306ac7' },
                { min: -20, max: -18, color: '#3b7edb' },
                { min: -20, max: -18, color: '#3b7edb' },
                { min: -18, max: -16, color: '#4e8add' },
                { min: -16, max: -14, color: '#6196e0' },
                { min: -14, max: -12, color: '#74a3e2' },
                { min: -12, max: -10, color: '#87afe5' },
                { min: -10, max: -8, color: '#9bbce8' },
                { min: -8, max: -6, color: '#9ac4dc' },
                { min: -6, max: -4, color: '#99cdd0' },
                { min: -4, max: -2, color: '#98d6c4' },
                { min: -2, max: 0, color: '#97e8ad' },
                { min: 0, max: 2, color: '#d7de7e' },
                { min: 2, max: 4, color: '#eadb70' },
                { min: 4, max: 6, color: '#f4d963' },
                { min: 6, max: 8, color: '#facc4f' },
                { min: 8, max: 10, color: '#f7b42d' },
                { min: 10, max: 12, color: '#f29b00' },
                { min: 12, max: 14, color: '#f19303' },
                { min: 14, max: 16, color: '#f0840a' },
                { min: 16, max: 18, color: '#ef7511' },
                { min: 18, max: 20, color: '#ee6618' },
                { min: 20, max: 22, color: '#ee581f' },
                { min: 22, max: 24, color: '#e74b1a' },
                { min: 24, max: 26, color: '#e03f16' },
                { min: 26, max: 28, color: '#d93312' },
                { min: 28, max: 30, color: '#d0240e' },
                { min: 30, max: 32, color: '#c20003' },
                { min: 32, max: 34, color: '#b50109' },
                { min: 34, max: 35, color: '#a90210' },
                { min: 35, max: 37, color: '#8a0519' },
                { min: 37, max: 40, color: '#6f0015' },
                { min: 40, color: '#6f0015' }
            ];
        var context, contours, d3Path, cnavas;
        function CanvasLayerInit(_data) {
            var canvasLayer = new ol.layer.Image({
                source: new ol.source.ImageCanvas({
                    canvasFunction: (extent, resolution, pixelRatio, size, projection) => {
                        const [width, height] = size;
                        const [left, bottom, right, top] = extent;
                        const xscale = width / (right - left);
                        const yscale = height / (top - bottom);

                        //获取裁剪区域的宽高
                        const pxregion = getPxRegion(Regioncoord, xscale, yscale, top, left);
                        const cw = Math.ceil(pxregion.xmax - pxregion.xmin)
                        const ch = Math.ceil(pxregion.ymax - pxregion.ymin)
                      
                        const regionTopLeft = pointToCoord(pxregion.xmin, pxregion.ymin, xscale, yscale, top, left);
                        const _left = regionTopLeft[0];
                        const _top = regionTopLeft[1];
                        const regionBottomRight = pointToCoord(pxregion.xmax, pxregion.ymax, xscale, yscale, top, left);
                        const _right = regionBottomRight[0];
                        const _bottom = regionBottomRight[1];

                        // 这里是canvas分辨率
                        const _xscale = cw / (_right - _left);
                        const _yscale = ch / (_top - _bottom);
                        const pxdata = getPxData(_data, _xscale, _yscale, _top, _left);

                        let canvas = document.createElement('canvas');
                        canvas.width = Math.ceil(pxregion.xmax);
                        canvas.height = Math.ceil(pxregion.ymax);
                        canvas.style.display = 'block';
                        // 设置canvas透明度
                        canvas.getContext('2d').globalAlpha = 0.1;
                        context = canvas.getContext('2d');
                        contours = d3.contours().size([cw, ch]); //等高线绘图实例
                        d3Path = d3.geoPath(null, context);  //绘图笔
                        context.clearRect(0, 0, cw, ch);

                        // 绘制裁剪区
                        context.beginPath()
                        for (let i = 0; i < pxregion.data.length; i++) {
                            const rp = pxregion.data[i]
                            context.lineTo(rp.x, rp.y)
                        }
                        context.closePath()
                        context.stroke()
                        context.clip()
                        context.globalAlpha = 0.3; //设置透明度
                        context.lineWidth = 2; //线条宽度
                        // 偏移
                        context.translate(pxregion.xmin, pxregion.ymin);

                        //进行空间插值
                        const idwdata = olIDW(pxdata.data, cw, ch);
                        //绘图等高线图
                        contours
                            .thresholds(10)
                            (idwdata)
                            .forEach(fill);
                        return canvas;
                    },
                    projection: 'EPSG:4326',
                    ratio: 1
                }),
                zIndex:15
            })
            // 向map添加图层
            mainMap.addLayer(canvasLayer);
        }
        function getPxRegion(data, xscale, yscale, top, left) {
            var _data = [];
            var xmax = 0;
            var ymax = 0;
            var xmin = 9999;
            var ymin = 9999;
            for (var i = 0; i < data.length; i++) {
                var p = coordToPoint(data[i][0], data[i][1], xscale, yscale, top, left)
                if (i === 0) {
                    xmax = p[0];
                    ymax = p[1];
                }
                if (p[0] != 0 || p[1] != 0) {
                    xmax > p[0] ? null : xmax = p[0];
                    ymax > p[1] ? null : ymax = p[1];
                    xmin < p[0] ? null : xmin = p[0];
                    ymin < p[1] ? null : ymin = p[1];
                    _data.push({
                        x: p[0],
                        y: p[1]
                    })
                }
            }
            return {
                data: _data,
                xmax: xmax,
                ymax: ymax,
                xmin: xmin,
                ymin: ymin,
            }
        }
        function getPxData(data, xscale, yscale, top, left) {
            var _data = [];
            var xmax = 0;
            var ymax = 0;
            var xmin = 9999;
            var ymin = 9999;
            for (var i = 0; i < data.length; i++) {
                var p = coordToPoint(data[i].Lon, data[i].Lat, xscale, yscale, top, left)
                if (i === 0) {
                    xmax = p[0];
                    ymax = p[1];
                }
                if (p[0] != 0 || p[1] != 0) {
                    xmax > p[0] ? null : xmax = p[0];
                    ymax > p[1] ? null : ymax = p[1];
                    xmin < p[0] ? null : xmin = p[0];
                    ymin < p[1] ? null : ymin = p[1];
                    _data.push({
                        x: p[0],
                        y: p[1],
                        value: Number(data[i].TEM)
                    })
                }
            }
            // value一定要确保是number类型，否则在插值计算的会先string转number再计算，影响效率
            return {
                data: _data,
                xmax: xmax,
                ymax: ymax,
                xmin: xmin,
                ymin: ymin,
            }
        }
        function pointToCoord(px, py, xscale, yscale, top, left) {
            const x = (px / xscale) + left;
            const y = top - (py / yscale);
            return [x, y];
        }
        function coordToPoint(x, y, xscale, yscale, top, left) {
            // 目的是为了剔除范围外的数据
            let px = (x - left) * xscale;
            let py = (top - y) * yscale;
            px > 0 ? null : px = 0;
            py > 0 ? null : py = 0;
            return [px, py]
        }
        function olIDW(data, width, height) {
            var s = new Date().getTime();
            var d = data;
            //已有点初始二维数组
            var dlen = d.length;
            var matrixData = []

            /**
            * 插值矩阵数据,时间复杂度O(height*width*len)
            */
            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 = 0; k < dlen; 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++
                    }
                    if (sum1 != 0)
                        //matrixData[k1] = sum0 / sum1 - referenceValue;
                        matrixData[k1] = sum0 / sum1;
                    else
                        matrixData[k1] = 1;
                }
            }
            var e = new Date().getTime();
            return matrixData;
        }
        //绘图一条等高线
        function fill(geometry) {
            context.beginPath();
            d3Path(geometry);
            var acolor = getColor(colors, geometry.value);
            //等高线梯度着色
            context.fillStyle = acolor;
            context.fill();
            //等高线线条
            context.strokeStyle = acolor;
            context.stroke();
        }
        function colorRGB(color) {
            let sColor = color.toLowerCase();
            // 处理六位的颜色值
            const sColorChange = [];
            // 十六进制颜色值的正则表达式
            const reg = /^#([0-9a-fA-f]{3}|[0-9a-fA-f]{6})$/;
            // 如果是16进制颜色
            if (sColor && reg.test(sColor)) {
                if (sColor.length === 4) {
                    let sColorNew = '#';
                    for (let i = 1; i < 4; i += 1) {
                        sColorNew += sColor.slice(i, i + 1).concat(sColor.slice(i, i + 1));
                    }
                    sColor = sColorNew;
                }
                for (let i = 1; i < 7; i += 2) {
                    sColorChange.push(parseInt('0x' + sColor.slice(i, i + 2)));
                }
            }
            return sColorChange;
        }
        function getColor(_colors, value) {
            var len = _colors.length;
            for (var i = 0; i < len; i++) {
                if (value > _colors[i].min && value <= _colors[i].max) {
                    const color = colorRGB(_colors[i].color);
                    return d3.rgb(color[0], color[1], color[2]);
                }
            }
        }
    </script>
</body>
</html>
参考文章：
1.反距离权重IDW算法解析与示例：https://www.jianshu.com/p/af471b72510d
2.反距离加权法（Inverse Distance Weighted）插值：https://www.jianshu.com/p/b38c5e464d16
3.使用D3-contour基于idw插值算法前端实现等值线色斑图：https://mp.weixin.qq.com/s/6BECmmW7evLHd2v018ffsA?client=tim&ADUIN=741841403&ADSESSION=1573171691&ADTAG=CLIENT.QQ.5603_.0&ADPUBNO=26933
**本文已结束，感谢您的观看**
鼓励一下

往期推荐

小额赞助
本人提供免费与付费咨询服务，感谢您的支持！赞助请发邮件通知，方便公布您的善意！
**光 3.01 元
Sun 3.00 元
bibichuan 3.00 元
捐赠者名单 合计：9.01 元
微信公众号
广告位
诚心邀请广大金主爸爸洽谈合作
每日一省
isNaN 和 Number.isNaN 函数的区别？

1.函数 isNaN 接收参数后，会尝试将这个参数转换为数值，任何不能被转换为数值的的值都会返回 true，因此非数字值传入也会返回 true ，会影响 NaN 的判断。

2.函数 Number.isNaN 会首先判断传入参数是否为数字，如果是数字再继续判断是否为 NaN ，不会进行数据类型的转换，这种方法对于 NaN 的判断更为准确。

2021-04-07 链接地址
每日二省
为什么0.1+0.2 ! == 0.3，如何让其相等?

一个直接的解决方法就是设置一个误差范围，通常称为“机器精度”。对JavaScript来说，这个值通常为2-52，在ES6中，提供了Number.EPSILON属性，而它的值就是2-52，只要判断0.1+0.2-0.3是否小于Number.EPSILON，如果小于，就可以判断为0.1+0.2 ===0.3。

2021-04-07 链接地址
每日三省
== 操作符的强制类型转换规则？

1.首先会判断两者类型是否**相同，**相同的话就比较两者的大小。

2.类型不相同的话，就会进行类型转换。

3.会先判断是否在对比 null 和 undefined，是的话就会返回 true。

4.判断两者类型是否为 string 和 number，是的话就会将字符串转换为 number。

5.判断其中一方是否为 boolean，是的话就会把 boolean 转为 number 再进行判断。

6.判断其中一方是否为 object 且另一方为 string、number 或者 symbol，是的话就会把 object 转为原始类型再进行判断。

2021-04-07 链接地址
每日英语
Happiness is time precipitation, smile is the lonely sad.
幸福是年华的沉淀，微笑是寂寞的悲伤。
2016 ——— 美文