Turf处理等压线

news2024/11/23 16:40:57

        Turf是一个用于空间分析的JavaScript库。它包括传统的空间操作、用于创建GeoJSON数据的辅助函数以及数据分类和统计工具。Turf可以作为客户端插件添加到您的网站,也可以使用Node.js在服务器端运行Turf。

        Turf是一个模块化的js类库,所有的模块都是在packages下面,每个模块下面都包含一个index.ts的源码文件,这个源码是使用typescript编写.

       

        对于上面这个实例,本章将通过源码的方式进行解析,首先调用的是turf.pointGrid这个函数,

这个模块中的index.ts文件,代码中定义了pointGrid函数,这个函数接受如下几个参数

// 基于box边界,通过细分的方式插入点
function pointGrid<P = GeoJsonProperties>(
  bbox: BBox,      // 东南西北作为边界
  cellSide: number,// 细分时的步长
  options: {        
    units?: Units; // 单位
    mask?: Feature<Polygon | MultiPolygon>; // 多边形掩码,插值点在这个范围内才算
    properties?: P;    // 属性
  } = {}
): FeatureCollection<Point, P> {
 // Default parameters
  // 默认使用千米
  if (options.mask && !options.units) options.units = "kilometers";

  // Containers
  // 结果
  var results = [];

  // 东南西北四个边界(经纬度)
  var west = bbox[0];
  var south = bbox[1];
  var east = bbox[2];
  var north = bbox[3];

  // 计算两点之间的球面距离

  // 一个单元占整个距离的百分比
  var xFraction = cellSide / distance([west, south], [east, south], options);
  // 一个单元的角度制下的距离
  var cellWidth = xFraction * (east - west);

  // 一个单元占整个距离的百分比
  var yFraction = cellSide / distance([west, south], [west, north], options);
  // 一个单元的角度制下的距离
  var cellHeight = yFraction * (north - south);

  // 宽度(角度)
  var bboxWidth = east - west;
  // 高度(角度)
  var bboxHeight = north - south;

  // 列数
  var columns = Math.floor(bboxWidth / cellWidth);
  // 行数
  var rows = Math.floor(bboxHeight / cellHeight);

  // adjust origin of the grid
  // 多出来的一部分的中间
  var deltaX = (bboxWidth - columns * cellWidth) / 2;
  var deltaY = (bboxHeight - rows * cellHeight) / 2;

  // 当前x
  var currentX = west + deltaX;

  // 计算插值点
  while (currentX <= east) {
    // 当前y
    var currentY = south + deltaY;
    while (currentY <= north) {
      // 封装要素点
      var cellPt = point([currentX, currentY], options.properties);
      if (options.mask) {
        // 点是否在掩码范围被内,只有在这个范围内才添加
        if (within(cellPt, options.mask)) results.push(cellPt);
      } else {
        // 添加到结果中
        results.push(cellPt);
      }
      currentY += cellHeight;
    }
    currentX += cellWidth;
  }

  // 封装成要素集合
  return featureCollection(results);
}

 计算distance的类是turf-distance文件夹下的index.ts

// 计算距离
function distance(
  from: Coord | Point,
  to: Coord | Point,
  options: {
    units?: Units;
  } = {}
) {
  // 获取坐标数组
  var coordinates1 = getCoord(from);
  // 获取坐标数组
  var coordinates2 = getCoord(to);
  // 计算维度
  var dLat = degreesToRadians(coordinates2[1] - coordinates1[1]);
  // 计算精度
  var dLon = degreesToRadians(coordinates2[0] - coordinates1[0]);
  // 角度转换为弧度
  var lat1 = degreesToRadians(coordinates1[1]);
  var lat2 = degreesToRadians(coordinates2[1]);

  // 计算两点之间的角度
  var a =
    Math.pow(Math.sin(dLat / 2), 2) +
    Math.pow(Math.sin(dLon / 2), 2) * Math.cos(lat1) * Math.cos(lat2);

  // 弧长公式Rθ
  return radiansToLength(
    2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)),
    options.units
  );
}

  计算等值线的类是turf-isolines文件夹下的index.ts

// 生成等值线
function isolines(
  pointGrid: FeatureCollection<Point>,// 网格点
  breaks: number[],
  options?: {
    zProperty?: string;
    commonProperties?: GeoJsonProperties;
    breaksProperties?: GeoJsonProperties[];
  }
) {
  // Optional parameters
  options = options || {};
  if (!isObject(options)) throw new Error("options is invalid");
  const zProperty = options.zProperty || "elevation";
  const commonProperties = options.commonProperties || {};
  const breaksProperties = options.breaksProperties || [];

  // Input validation
  collectionOf(pointGrid, "Point", "Input must contain Points");
  if (!breaks) throw new Error("breaks is required");
  if (!Array.isArray(breaks)) throw new Error("breaks must be an Array");
  if (!isObject(commonProperties))
    throw new Error("commonProperties must be an Object");
  if (!Array.isArray(breaksProperties))
    throw new Error("breaksProperties must be an Array");

  // Isoline methods
  // 等值线方法

  // 将网格点,按照z值构建一个只包含z值的矩阵
  const matrix = gridToMatrix(pointGrid, { zProperty: zProperty, flip: true });
  // 创建等值线
  const createdIsoLines = createIsoLines(
    matrix,
    breaks,
    zProperty,
    commonProperties,
    breaksProperties
  );
  // 重新缩放等值线
  const scaledIsolines = rescaleIsolines(createdIsoLines, matrix, pointGrid);

  // 封装成要素
  return featureCollection(scaledIsolines);
}

将所有的点构建一个z值的矩阵

export default function gridToMatrix(grid, options) {
  // Optional parameters
  options = options || {};
  if (!isObject(options)) throw new Error("options is invalid");
  var zProperty = options.zProperty || "elevation";
  var flip = options.flip;
  var flags = options.flags;

  // validation
  collectionOf(grid, "Point", "input must contain Points");

  // 按照经纬度排序
  var pointsMatrix = sortPointsByLatLng(grid, flip);

  var matrix = [];
  // create property matrix from sorted points
  // looping order matters here
  // 从排序点创建属性矩阵循环此处的顺序问题
  for (var r = 0; r < pointsMatrix.length; r++) {
    var pointRow = pointsMatrix[r];
    var row = [];
    for (var c = 0; c < pointRow.length; c++) {
      var point = pointRow[c];
      // Check if zProperty exist
      if (point.properties[zProperty]) row.push(point.properties[zProperty]);
      else row.push(0);
      // add flags
      if (flags === true) point.properties.matrixPosition = [r, c];
    }
    matrix.push(row);
  }

  return matrix;
}

 

// 创建等值线
function createIsoLines(
  matrix: number[][],
  breaks: number[],
  zProperty: string,
  commonProperties: GeoJsonProperties,
  breaksProperties: GeoJsonProperties[]
): Feature<MultiLineString>[] {
  const results = [];
  // 遍历中断点
  for (let i = 1; i < breaks.length; i++) {
    // 作为阈值
    const threshold = +breaks[i]; // make sure it's a number

    const properties = { ...commonProperties, ...breaksProperties[i] };
    properties[zProperty] = threshold;
    // 应该是按照矩阵中的四个点计算是否可以计算出这个阈值的插值点
    const isoline = multiLineString(isoContours(matrix, threshold), properties);

    results.push(isoline);
  }
  return results;
}
// 等压线
export default function isoContours(data, threshold, options) {
  /* process options */
  options = options ? options : {};

  var optionKeys = Object.keys(defaultSettings);

  for (var i = 0; i < optionKeys.length; i++) {
    var key = optionKeys[i];
    var val = options[key];
    val =
      typeof val !== "undefined" && val !== null ? val : defaultSettings[key];

    settings[key] = val;
  }

  if (settings.verbose)
    console.log(
      "MarchingSquaresJS-isoContours: computing isocontour for " + threshold
    );

    // 计算等压线网格
  var ret = contourGrid2Paths(computeContourGrid(data, threshold));

  if (typeof settings.successCallback === "function")
    settings.successCallback(ret);

  return ret;
}
/* assume that x1 == 1 &&  x0 == 0 */
// 线性插值
function interpolateX(y, y0, y1) {
  return (y - y0) / (y1 - y0);
}

/* compute the isocontour 4-bit grid */
// 计算等压线4为的网格
function computeContourGrid(data, threshold) {
  // 行数、列数
  var rows = data.length - 1;
  var cols = data[0].length - 1;
  // 
  var ContourGrid = { rows: rows, cols: cols, cells: [] };

  // 遍历行数
  for (var j = 0; j < rows; ++j) {
    ContourGrid.cells[j] = [];
    // 遍历列数
    for (var i = 0; i < cols; ++i) {
      /* compose the 4-bit corner representation */
      var cval = 0;

      //取出相邻的四个点
      var tl = data[j + 1][i];
      var tr = data[j + 1][i + 1];
      var br = data[j][i + 1];
      var bl = data[j][i];

      if (isNaN(tl) || isNaN(tr) || isNaN(br) || isNaN(bl)) {
        continue;
      }

      // 检索阈值是否在四个高度内部,还是外部,四个值都高于阈值8421=15,或者四个值都低于阈值在外部0
      cval |= tl >= threshold ? 8 : 0;
      cval |= tr >= threshold ? 4 : 0;
      cval |= br >= threshold ? 2 : 0;
      cval |= bl >= threshold ? 1 : 0;

      /* resolve ambiguity for cval == 5 || 10 via averaging */
      var flipped = false;
      // 在对角线上
      if (cval === 5 || cval === 10) {
        // 计算平均值
        var average = (tl + tr + br + bl) / 4;
        // 两个对角线上的中点值偏离比较远,
        if (cval === 5 && average < threshold) {
          cval = 10; // 取另外一个对角线来确定
          flipped = true;
        } else if (cval === 10 && average < threshold) {
          cval = 5;
          flipped = true;
        }
      }

      /* add cell to ContourGrid if it contains edges */
      // 如果单元格包含边,则将其添加到ContourGrid

      // 阈值在四个点内部
      if (cval !== 0 && cval !== 15) {

        var top, bottom, left, right;
        top = bottom = left = right = 0.5;
        /* interpolate edges of cell */
        if (cval === 1) {  // 一个点在上,其他三个点都在下面
          // 插值左边
          left = 1 - interpolateX(threshold, tl, bl);
          // 插值下边
          bottom = 1 - interpolateX(threshold, br, bl);
        } else if (cval === 2) {
          bottom = interpolateX(threshold, bl, br);
          right = 1 - interpolateX(threshold, tr, br);
        } else if (cval === 3) {
          left = 1 - interpolateX(threshold, tl, bl);
          right = 1 - interpolateX(threshold, tr, br);
        } else if (cval === 4) {
          top = interpolateX(threshold, tl, tr);
          right = interpolateX(threshold, br, tr);
        } else if (cval === 5) {
          top = interpolateX(threshold, tl, tr);
          right = interpolateX(threshold, br, tr);
          bottom = 1 - interpolateX(threshold, br, bl);
          left = 1 - interpolateX(threshold, tl, bl);
        } else if (cval === 6) {
          bottom = interpolateX(threshold, bl, br);
          top = interpolateX(threshold, tl, tr);
        } else if (cval === 7) {
          left = 1 - interpolateX(threshold, tl, bl);
          top = interpolateX(threshold, tl, tr);
        } else if (cval === 8) {
          left = interpolateX(threshold, bl, tl);
          top = 1 - interpolateX(threshold, tr, tl);
        } else if (cval === 9) {
          bottom = 1 - interpolateX(threshold, br, bl);
          top = 1 - interpolateX(threshold, tr, tl);
        } else if (cval === 10) {
          top = 1 - interpolateX(threshold, tr, tl);
          right = 1 - interpolateX(threshold, tr, br);
          bottom = interpolateX(threshold, bl, br);
          left = interpolateX(threshold, bl, tl);
        } else if (cval === 11) {
          top = 1 - interpolateX(threshold, tr, tl);
          right = 1 - interpolateX(threshold, tr, br);
        } else if (cval === 12) {
          left = interpolateX(threshold, bl, tl);
          right = interpolateX(threshold, br, tr);
        } else if (cval === 13) {
          bottom = 1 - interpolateX(threshold, br, bl);
          right = interpolateX(threshold, br, tr);
        } else if (cval === 14) {
          left = interpolateX(threshold, bl, tl);
          bottom = interpolateX(threshold, bl, br);
        } else {
          console.log(
            "MarchingSquaresJS-isoContours: Illegal cval detected: " + cval
          );
        }
        ContourGrid.cells[j][i] = {
          cval: cval,   // 代表的是top、right、bottom、left中哪两个值有效果
          flipped: flipped,
          top: top,
          right: right,
          bottom: bottom,
          left: left,
        };
      }
    }
  }

  return ContourGrid;
}
// 计算等值点组成的路径
function contourGrid2Paths(grid) {
  var paths = [];
  var path_idx = 0;
  var epsilon = 1e-7;

  grid.cells.forEach(function (g, j) {
    g.forEach(function (gg, i) {
      if (typeof gg !== "undefined" && !isSaddle(gg) && !isTrivial(gg)) {
        var p = tracePath(grid.cells, j, i);
        var merged = false;
        /* we may try to merge paths at this point */
        if (p.info === "mergeable") {
          /*
            search backwards through the path array to find an entry
            that starts with where the current path ends...
          */
          var x = p.path[p.path.length - 1][0],
            y = p.path[p.path.length - 1][1];

          for (var k = path_idx - 1; k >= 0; k--) {
            if (
              Math.abs(paths[k][0][0] - x) <= epsilon &&
              Math.abs(paths[k][0][1] - y) <= epsilon
            ) {
              for (var l = p.path.length - 2; l >= 0; --l) {
                paths[k].unshift(p.path[l]);
              }
              merged = true;
              break;
            }
          }
        }
        if (!merged) paths[path_idx++] = p.path;
      }
    });
  });

  return paths;
}

总结:

具体的做法是按照给定的采样点,插值(例如反距离权重)成矩阵,将高度值构建成矩阵,将矩阵中的四个相邻的值取出来,按照给定的阈值进行线性插值定位,计算z值对应的经纬度,这就构成一个插值后的xyz值,将这些值按照路径拼起来就是等直线了

参考:

等值线 | Turf.js中文网

GitHub - Turfjs/turf: A modular geospatial engine written in JavaScript

等值线 | Turf.js中文网

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1051268.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

【GDB】用 python 扩展 gdb

用 python 扩展 GDB .gdbinit 文件中实现自定义命令 mv 代码如下 define mvif $argc 2delete $arg0# 注意新创建的断点编号和被删除断点的编号不同break $arg1elseprint "输入参数数目不对&#xff0c;help mv 以获得用法"end end# (gdb) help mv 会输出以下帮助文…

搭建自己的搜索引擎之四

一、前言 搭建自己的搜索引擎之三 介绍了通过HTTP RESTful 对ES进行增删改查&#xff0c;这一般手工运维ES时使用&#xff0c;程序代码中最好还是使用Java API去操作ES会更容易维护&#xff0c;但ES API竟然贼多&#xff0c;本篇介绍一下 四种 API及其简单使用。 注&#xff…

深入理解二叉树:结构、遍历和实现

文章目录 &#x1f34b;引言&#x1f34b;什么是二叉树&#xff1f;&#x1f34b;二叉树的基本性质&#x1f34b;二叉树的遍历&#x1f34b;二叉树的实现&#x1f34b;结语 &#x1f34b;引言 在计算机科学中&#xff0c;二叉树是一种重要的数据结构&#xff0c;广泛应用于各种…

【密评】商用密码应用安全性评估从业人员考核题库(二)

商用密码应用安全性评估从业人员考核题库&#xff08;二&#xff09; 国密局给的参考题库5000道只是基础题&#xff0c;后续更新完5000还会继续更其他高质量题库&#xff0c;持续学习&#xff0c;共同进步。 251 多项选择题 根据《密码法》&#xff0c;核心密码、普通密码安全…

Linux常用命令(一)

目录 一、列出目录内容&#xff08;ls&#xff09; 二、切换目录&#xff08;cd&#xff09; 三、显示当前目录路径&#xff08;pwd&#xff09; 四、以树状结构显示目录内容&#xff08;tree&#xff09; 五、创建新目录&#xff08;mkdir&#xff09; 六、复制文件或目…

windows11 如何关闭 vbs

在Windows 11中&#xff0c;VBS是一种虚拟化安全功能&#xff0c;它可以防止恶意软件通过沙箱环境运行。 如果您想关闭VBS功能&#xff0c;方法如下&#xff1a; 点击底部开始菜单 在上方搜索 cmd &#xff0c;并点击以管理员身份运行 打开控制台后&#xff0c;在控制台输入…

文档图像处理:大模型的突破与新探索

前言 随着数字化时代的到来&#xff0c;文档图像处理技术在各行各业扮演着越来越重要的角色。在2023第十二届中国智能产业高峰论坛&#xff08;CIIS 2023&#xff09;的专题论坛上&#xff0c;合合信息智能技术平台事业部副总经理、高级工程师丁凯博士分享了当前文档图像处理面…

wallis匀色算法、直方图匹配、颜色转移方法比较

算法原理 这三种方法应该是比较基础的匀色处理算法 三个算法的原理比较简单&#xff0c;具体原理大家可以自己百度 &#xff08;1&#xff09;wallis匀色原理主要在于利用Wallis滤波器使原始图像的均值和标准差与参考影像相当&#xff0c;从而使原始影像和参考影像具有相近的色…

Oracle的递归公共表表达式

查询节点id为2的所有子节点的数据&#xff0c;包括向下级联 WITH T1 (id, parent_id, data) AS (SELECT id, parent_id, dataFROM nodesWHERE id 2UNION ALLSELECT t.id, t.parent_id, t.dataFROM nodes tJOIN T1 n ON t.parent_id n.id ) SELECT * FROM T1; --建表语句 C…

今天出门竟然忘了带套

今天是没有抢到票的打工人节前的最后一天&#xff0c;7点醒来&#xff0c;磨磨蹭蹭&#xff0c;解决完个人问题&#xff0c;7.35才出门&#xff0c;正常来说最晚7.30出门&#xff0c;骑上哈啰、挤上地铁才能保证打上卡。 说出来不怕各位同行笑话&#xff0c;谁能想到一个高速发…

打卡新“姿势”,多种打卡方式并行

打卡工具 路径 拓展 >> 工具 功能简介 在打卡工具 “班次管理”中&#xff0c;支持多种打卡方式。可同时选择「地点打卡」和「智能安全帽打卡」两种方式进行打卡。 注&#xff1a; 「地点打卡」可设置考勤地点&#xff1b; 「智能安全帽打卡」可设置电子围栏范围。…

排序篇(一)----插入排序

1.直接插入排序 插入排序的思想: 把待排序的记录按其关键码值的大小逐个插入到一个已经排好序的有序序列中&#xff0c;直到所有的记录插入完为止&#xff0c;得到一个新的有序序列 。 你可以想像成打牌一样,比如说斗地主,一张一张的摸牌,然后把手上的这些牌变成手续的排列.…

【教学类-38】A4红纸-国旗灯笼(庆祝中华人民共和国成立74周年)

作品展示&#xff1a; 背景需求&#xff1a; 从教十余年&#xff0c;我在每年国庆都带领中大班孩子们制作与“国旗相关”国庆庆祝物品——国旗、礼盒 一、国旗&#xff08;吸管、A4红纸、黄纸打印五角星&#xff09; 二、铅画纸手提袋&#xff08;8K铅画纸、A4红纸、黄色打印…

Windows的批处理——获取系统时间、生成当天日期日志

Windows批处理基础https://coffeemilk.blog.csdn.net/article/details/132118351 一、Windows批处理的日期时间 在我们进行软件开发的过程中&#xff0c;有时候会使用到一些批处理命令&#xff0c;其中就涉及到获取系统日期、时间来进行一些逻辑的判断处理&#xff1b;那么我们…

Ubuntu 部署 Seata1.7.1

一、环境说明 IP操作系统程序备注10.0.61.22ubuntu20.04PostgreSQL-14.11已提前部署10.0.61.21ubuntu20.04Nacos-2.1.0已提前部署10.0.61.22ubuntu20.04seata-server-1.7.1本文将要部署 二、部署 1. 下载 wget https://github.com/seata/seata/releases/download/v1.7.1/se…

VUE2项目:尚品汇VUE-CLI脚手架初始化项目以及路由组件分析(一)

标题 环境VUE2目录publicassetscomponentsmain.jsbabel.config.jspackage.jsonvue.config.js 项目路由分析Header与Footer非路由组件完成Header示例 路由组件的搭建声明式导航编程式导航 Footer组件的显示与隐藏路由传递参数重写push和replace三级联动组件拆分附件 环境 前提要…

Scala第四章节

Scala第四章节 scala总目录 章节目标 掌握分支结构的格式和用法掌握for循环和while循环的格式和用法掌握控制跳转语句的用法掌握循环案例理解do.while循环的格式和用法 1. 流程控制结构 1.1 概述 在实际开发中, 我们要编写成千上万行代码, 代码的顺序不同, 执行结果肯定也…

GD32工程创建

1.创建空工程 在任意路径下创建空的test文件夹。打开keil5空工程创建空工程 选择对应的芯片型号&#xff1a; 然后把空工程保存到test文件夹下。会自动生成如下文件。 2. 添加组 下载GD32F10X的固件库&#xff1a;在百度里搜索GD32进入官网。 下载下来对应的文件如下&#xff…

问题记录 springboot 事务方法中使用this调用其它方法

原因: 因为代理对象中调用了原始对象的toString()方法,所以两个不同的对象打印出的引用是相同的

HTML详细基础(三)表单控件

本帖介绍web开发中非常核心的标签——表格标签。 在日常我们使用到的各种需要输入用户信息的场景——如下图&#xff0c;均是通过表格标签table创造出来的&#xff1a; 目录 一.表格标签 二.表格属性 三.合并单元格 四.无序列表 五.有序列表 六.自定义标签 七.表单域 …