5.gis空间插值之前端turf+kriging

最新推荐文章于 2024-04-30 22:33:50 发布

紫雪giser

最新推荐文章于 2024-04-30 22:33:50 发布

阅读量317

点赞数 1

分类专栏： gis开发学习文章标签：学习方法 arcgis 开源

本文链接：https://blog.csdn.net/yaojindou0/article/details/132038280

版权

gis开发学习专栏收录该内容

32 篇文章 8 订阅

订阅专栏

前言

本文主要讲述通过前端来实现通过空间插值生成等值线面、泰森多边形、TIN三角网,除此也可使用arcgis发布gp模型服务，再通过arcgis api调用

参数准备

import kriging from '@sakitam-gis/kriging'
import * as turf from '@turf/turf'

let params = {
    krigingModel: 'exponential',//model还可选'gaussian'高斯,'spherical'球状,exponential指数
    krigingSigma2: 0,
    krigingAlpha: 100,
    poygon: {//等值面
        id: {
           
            colors: ["#00a6ff", "#00ebff", "#00ffb5", "#00ff94", "#00ff08", "#10ff00", "#31fe00", "#52fe00","#9bff00", "#94ff00", "#b5ff00", "#deff00", "#ffff00", "#ffe700", "#ffc300","#ffa200", "#ff8200", "#ff6100", "#ff4100", "#ff2c00", "#ff2000", "#ff0000", "#790024", "#790024"],
            breaks: [-20, -16, -12, -8, -4, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 33, 36, 40],
            values: {
                "-20--16": 1, "-16--12": 2, "-12--8": 3, "-8-0": 4, "0-2": 5, "2-4": 6, "4-6": 7,
                "6-8": 8, "8-10": 9, "10-12": 10, "12-14": 11, "14-16": 12, "16-18": 13, "18-20": 14,
                "20-22": 15, "22-24": 16, "24-26": 17, "26-28": 18, "28-30": 19, "30-33": 20, "33-36": 21, "36-40": 22
            }
        }

    },
     line: {//等值线
        sx_dq_qxzqyid: {//气压
            colors: ["#96c69a", "#45c681", "#04a95a", "#04d209", "#070fde"],
            values: { "800": 0, "810": 0, "820": 1, "840": 1, "860": 1, "880": 2, "900": 2, "920": 3, "940": 3, "950": 4, "960": 4, "970": 4, "980": 4, "1200": 4 },
            breaks: [800, 810, 820, 840, 860, 880, 900, 920, 940, 950, 960, 970, 980, 1200]
        }
    }
 }

克里金插值


// geom 裁剪范围坐标 items点集[{longitude:108,latitude:34,value}] type线面类型 filed字段名
export const drawKriging = (geom, items, type, filed) => {
   let values = [], lngs = [], lats = [];
   items.forEach(item => {
       values.push(item[filed]);
       lngs.push(item.longitude);
       lats.push(item.latitude);
   });

   if (values.length > 3) {
       //kriging插值生成格网点集
       let variogram = kriging.train(values, lngs, lats, params.krigingModel, params.krigingSigma2, params.krigingAlpha);
       let polygons = geom.coordinates[0];
       let extent = geom.extent;

       let grid = kriging.grid(polygons, variogram, (extent[2] - extent[0]) / 100);
       let fc = gridFeatureCollection(grid);
       var collection = turf.featureCollection(fc);
       var features;
       let breaks = params[type][layerid].breaks;
       if (type === "poygon") {
           // 等值面生成
           features = turf.isobands(collection, breaks, { zProperty: 'value' });
           const sortArea = (a, b) => {
               return turf.area(b) - turf.area(a);
           };
           //按照面积对图层进行排序，规避turf的一个bug
           features.features.sort(sortArea)
       } else if (type === "line") {
           // 等值线生成
           features = turf.isolines(collection, breaks, { zProperty: 'value' });
       }

      
   } else {
       alert("有效样点个数不足，无法插值");
   }
}

反距离权重

// geom 裁剪范围坐标 items点集[{longitude:108,latitude:34,value}] type线面类型 filed字段名
export const drawFjlqz = (geom, items, type, filed) => {
    let bbs = [];
    items.forEach(item => {
        bbs.push(turf.point([item.longitude, item.latitude], { value: item[filed] }));
    });
    if (bbs.length > 3) {
        var points = turf.featureCollection(bbs);
        // 反距离权重插值生成格网点集
          let polygons = geom.coordinates[0];
        let extent = geom.extent;
        var collection = turf.interpolate(points, (extent[2] - extent[0]) / 100, { gridType: 'points', property: 'value', units: 'degrees' });
        var features;
        let breaks = params[type][layerid].breaks;
        if (type === "poygon") {
            // 等值面生成
            features = turf.isobands(collection, breaks, { zProperty: 'value' });
            const sortArea = (a, b) => {
                return turf.area(b) - turf.area(a);
            };
            //按照面积对图层进行排序，规避turf的一个bug
            features.features.sort(sortArea)
        } else if (type === "line") {
            // 等值线生成
            features = turf.isolines(collection, breaks, { zProperty: 'value' });
        }
        let intersectfeatures = []
        let polygon = turf.polygon(coord);
        features.features.forEach(feature => {
            let intersection = turf.intersect(polygon, feature);
            if (intersection) {
                intersection.properties = feature.properties;
                intersectfeatures.push(intersection);
            }
        });
        let intersections = turf.featureCollection(intersectfeatures);
       
    } else {
        alert("有效样点个数不足，无法插值");
    }
}

泰森多边形

// geom裁剪范围坐标 items点集[{longitude:108,latitude:34,value}]  filed字段名
export const drawTsdbx = (geom, items, filed) => {
    let bbs = [];
    items.forEach(item => {
        bbs.push(turf.point([item.longitude, item.latitude], { value: item[filed] }));
    });
    if (bbs.length > 3) {
        var points = turf.featureCollection(bbs);
         let coord = geom.coordinates[0];
        let extent = geom.extent;
        var features = turf.voronoi(points, { bbox: extent });
        let intersectfeatures = []
        let polygon = turf.polygon(coord);
        features.features.forEach(feature => {
            let intersection = turf.intersect(polygon, feature);
            if (intersection) {
                intersection.properties = feature.properties;
                intersectfeatures.push(intersection);
            }
        });
        let intersections = turf.featureCollection(intersectfeatures);
    
    } else {
        alert("有效样点个数不足，无法计算");
    }
}

TIN三角网

// geom 裁剪范围坐标 items点集[{longitude:108,latitude:34,value}] filed字段名
export const drawTin = (geom, items, filed) => {
    let bbs = [];
    items.forEach(item => {
        bbs.push(turf.point([item.longitude, item.latitude], { value: item[filed] }));
    });
    if (bbs.length > 3) {
        var points = turf.featureCollection(bbs);
         let coord = geom.coordinates[0];
      
        var features = turf.tin(points, 'value');
        let intersectfeatures = []
        let polygon = turf.polygon(coord);
        features.features.forEach(feature => {
            let intersection = turf.intersect(polygon, feature);
            if (intersection) {
                intersection.properties = feature.properties;
                intersectfeatures.push(intersection);
            }
        });
        let intersections = turf.featureCollection(intersectfeatures);
      
    } else {
        alert("有效样点个数不足，无法计算");
    }
}

总结

IDW 法概念很简单，每个内插点 (网格) 的值与邻近样本点的关系是距离，距离越远关系越小。所以取值点与样本点间的距离为权重进行加权平均，离内插点越近的样本点赋予的权重越大。

TIN 是不规则三角网，在电脑视觉领域或是GIS都很常见，跟前面 IDW 法是内插在规则的形状中有所不同，TIN 组成的不规则三角网是将样本点连成连续的三角网，而在众多产生三角网的算法中，Delaunay 三角化是公认最佳解：
Delaunay 三网化：数据中任三点取其外置圆，若此圆内没有包含任何其它点，则这三角形加入三角网中。这样的目的是让三角形都能越接近正三角形，狭长得三角形出现机会越低，因为三角形三边长若越接近，外置圆越小。
另外，GIS 领域有很多其他更多内插方法，例如 Bilinear, Kriging, nearest neighbor 等等