Kriging插值并生成Topojson

标签： Gis 分类： Gis 创建时间：2019-11-23 02:18:10 更新时间：2024-11-15 10:49:43

最开始的想法是从Free_GIS，公瑾的想法中产生的。原理使用开源的Kriging.js生成插值数据，然后再使用d3.contour进行等值线绘制，将坐标从笛卡尔坐标转换到地理坐标，最后生成GeoJSON。有些数据效果还是不错的。

这个算法主要有两个问题，第一就是在进行坐标转换的时候，他仅仅考虑了将像素点以图像的中心点旋转90度，没有考虑到长宽比不是一的情况，没有进行缩放。第二个问题，就是部分数据生成的数据绘制不全。像下面这组数据：

1
2
3
4
5
6
7
8
9
10
11
12
13
dataset={"type":"FeatureCollection","features":[{"type":"Feature","geometry":
        {"type":"Point","coordinates":[119.84338055,28.39903371]},"properties":{"id":40,"level":12},"id":"hlgy_station.1"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84262228,28.39976446]},"properties":{"id":57,"level":28},"id":"hlgy_station.2"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84365151,28.39991672]},"properties":{"id":58,"level":60},"id":"hlgy_station.3"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84268319,28.39863171]},"properties":{"id":59,"level":58},"id":"hlgy_station.4"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84309731,28.39874133]},"properties":{"id":80,"level":47},"id":"hlgy_station.5"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84231169,28.39945996]},"properties":{"id":81,"level":75},"id":"hlgy_station.6"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84238477,28.39940515]},"properties":{"id":85,"level":9},"id":"hlgy_station.7"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84246394,28.39934425]},"properties":{"id":92,"level":46},"id":"hlgy_station.8"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84311558,28.40037956]},"properties":{"id":82,"level":18},"id":"hlgy_station.9"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84299378,28.3994965]},"properties":{"id":83,"level":30},"id":"hlgy_station.10"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84306686,28.39944169]},"properties":{"id":84,"level":19},"id":"hlgy_station.11"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84393165,28.39923463]},"properties":{"id":94,"level":55},"id":"hlgy_station.12"}]}

就会出现数据绘制不全的情况。

无奈大神比较忙，没空搭理我。所以只能自己动手解决这些毛病了。

使用原始的kriging.js生成栅格图像的过程，还是比较顺利的。

问题主要出在将栅格图像进行等值线话的过程。关于坐标旋转的问题，我觉得需要使用空间变换中的仿射变换比较好。

参考文章：
1.FreeGIS/kriging-contour
2.Javascript图像处理——仿射变换
3.OpenLayers实现克里金插值渲染图-Kriging
4.Openlayers利用kriging.js实现纯前端插值
5.kriging.js

后来我尝试了使用turf进行生成等线的方法。经过谷歌，发现了一个好的网站，就是参考文章5，这里完全实现了通过kriging生成插值，并通过turf生成等值线的代码，然后转换成geojson格式加载到地图上的，我下下来，虽然和原始的kriging生成的栅格有些差距，也有些许的不如意，但是总体还是满意的。

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
var restify = require('restify');
var turf=require("@turf/turf");
var kriging=require("./lib/kriging");
const topojson=require("topojson");
var turf=require("@turf/turf");

const server = restify.createServer({
    name: 'PhGIShlgyNode',
    version: '1.0.0'
  });
  
server.use(restify.plugins.acceptParser(server.acceptable));
server.use(restify.plugins.queryParser());
server.use(restify.plugins.bodyParser());

server.use(function(req,res,next){
    res.header("Access-Control-Allow-Origin", "*");
    res.header("Access-Control-Allow-Headers", "Content-Type,Content-Length, Authorization, Accept,X-Requested-With");
    res.header("Access-Control-Allow-Methods","PUT,POST,GET,DELETE,OPTIONS");
    next();
})

// 权重生成等值面数组
/**
 * freegis.kriging_contour(dataset,weight_field,kriging_params,weight_breaks);
 * 
 * dataset：geojson格式的featureclass数据集，feature是图形是点
 * weight_field：绑定权重字段名称
 * kriging_params：克里金插值参数
 * weight_breaks：权重生成等值面分级数组
 * 
 */
let weight_breaks=[];
// 设置最大最小值，以及颜色变化的步长
let colorSize=20;
let maxValue=5000;
let minValue=50;
let step=(maxValue-minValue)/colorSize-2;

weight_breaks.push(minValue);
for(let i=1;i<colorSize-1;i++){
    weight_breaks.push(minValue+i*step);
}
weight_breaks.push(maxValue);

server.post('/api/kriging', function (req, res, next) {

    let dataset={"type":"FeatureCollection","features":[{"type":"Feature","geometry":
        {"type":"Point","coordinates":[119.84338055,28.39903371]},"properties":{"id":40,"value":12},"id":"hlgy_station.1"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84262228,28.39976446]},"properties":{"id":57,"value":28},"id":"hlgy_station.2"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84365151,28.39991672]},"properties":{"id":58,"value":60},"id":"hlgy_station.3"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84268319,28.39863171]},"properties":{"id":59,"value":58},"id":"hlgy_station.4"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84309731,28.39874133]},"properties":{"id":80,"value":47},"id":"hlgy_station.5"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84231169,28.39945996]},"properties":{"id":81,"value":75},"id":"hlgy_station.6"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84238477,28.39940515]},"properties":{"id":85,"value":9},"id":"hlgy_station.7"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84246394,28.39934425]},"properties":{"id":92,"value":46},"id":"hlgy_station.8"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84311558,28.40037956]},"properties":{"id":82,"value":18},"id":"hlgy_station.9"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84299378,28.3994965]},"properties":{"id":83,"value":30},"id":"hlgy_station.10"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84306686,28.39944169]},"properties":{"id":84,"value":19},"id":"hlgy_station.11"},
        {"type":"Feature","geometry":{"type":"Point","coordinates":[119.84393165,28.39923463]},"properties":{"id":94,"value":55},"id":"hlgy_station.12"}]}

    let values = [], lngs = [], lats = [];
    dataset.features.forEach(feature => {
        values.push(feature.properties.value);
        lngs.push(feature.geometry.coordinates[0]);
        lats.push(feature.geometry.coordinates[1]);
    });
    let variogram = kriging.train(
        values,
        lngs,
        lats,
        "exponential",
        0,
        100
    );
    let extent=[119.83800000000001, 28.395, 119.848, 28.404999999999998];

    let polygons = [];
    polygons.push([
        [extent[0], extent[1]], [extent[0], extent[3]],
        [extent[2], extent[3]], [extent[2], extent[1]]
    ]);
    let grid = kriging.grid(polygons, variogram, (extent[2] - extent[0]) / 500);
    let fc = gridFeatureCollection(grid,
        [extent[0], extent[2]], [extent[1], extent[3]]);
    var collection = turf.featureCollection(fc);
    var breaks = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0];
    var isobands = turf.isobands(collection, breaks, { zProperty: 'value' });
    function sortArea(a,b)
    {
        return turf.area(b)-turf.area(a);
    }
    //按照面积对图层进行排序，规避turf的一个bug
    isobands.features.sort(sortArea);

    res.send(isobands);
    return next();
});

server.listen(8799, function () {
    console.log('%s listening at %s', server.name, server.url);
});
//利用网格计算点集
const gridFeatureCollection = function (grid, xlim, ylim) {
    var range =grid.zlim[1] - grid.zlim[0];
    var i, j, x, y, z;
    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] - grid.zlim[0]) / range;
            if (z < 0.0) z = 0.0;
            if (z > 1.0) z = 1.0;
            pointArray.push(turf.point([x, y], { value: z }));
        }
    return pointArray;
}

效果如下图：

当然还可以进行平滑处理，就像参考文章3中的一样，但我觉得可以有，但是没有必要了。

参考文章：
1.OL4结合turf.js实现等值线
2.浅谈WebGIS等值线模态前端展示生成绘制(算法经验)
3.web前端利用turf.js生成等值线、等值面
4.Turfjs
5.OpenLayers 5 使用turf.js渲染克里金插值计算的等值面 (别担心，我打开了，是个安全的网站)

更新
经过长时间探索,我将差值后生成的geojson进行了裁切，而且使用了turf.simplify进行了抽稀，最后讲geojson转化为了topojson，大大的减轻了传输数据，而且也没有非常明显的锯齿效果。

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
// 前面的代码省略了

//是否需要裁切，裁切范围也是一个geojson
if(clip_geom){
    isobands.features.forEach(feature=>{
        var coors=martinez.intersection(feature.geometry.coordinates,clip_geom.coordinates);
        feature.geometry.coordinates=coors;
        return feature;
    });
}

//对Geojson进行简化
var resultFeatures=isobands.features;
var fcout=resultFeatures.length;
for(let i=fcout-1;i>=0;i--){
    let feature=resultFeatures[i].geometry;
    if(feature.coordinates.length>0){
        let options = {tolerance: 0.00000098, highQuality: false};
        let simplified = turf.simplify(feature, options);
        resultFeatures[i].geometry=simplified;
    }
}
// 将GeoJSON转为TopoJSON
isobands = topojson.topology({foo: isobands});

// 后面的代码省略了

**本文已结束，感谢您的观看**
鼓励一下

往期推荐

小额赞助
本人提供免费与付费咨询服务，感谢您的支持！赞助请发邮件通知，方便公布您的善意！
**光 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 ——— 美文