【NetTopologySuite类库】常用功能整理(1)

介绍

简称NTS,是C#端的GIS相关的拓扑库。Github地址API文档

其他相关

关于OGC标准(开放地理空间信息联盟(Open Geospatial Consortium))

项目环境

  1. 本文以VS2019创建的.Net Framework4.7.2 控制台应用为例,对NTS的一些功能使用进行介绍。
  2. 右击项目的引用,打开NuGet程序包管理器,搜索NetTopologySuite,安装NetTopologySuite1.14.0版和NetTopologySuite.IO1.14.0.1版本。

如果要看源码,考虑到网速问题,可以去Gitee上(已从Github中迁移过来,访问地址),根据分支中的标签,选择与NuGet中相同版本的源码。
在这里插入图片描述

常用的API介绍

IO相关

读取wkt字符串

作用:WKTReaderWKT格式的字符串转化为Geometry的对象。

例子:

using NetTopologySuite.Geometries;
using NetTopologySuite.IO;

namespace NTSRecipes
{
    class Program
    {
        static void Main(string[] args)
        {
            WKTReader reader = new WKTReader();
            var point = reader.Read("Point(0 0)") as Point;
            var line = reader.Read("LineString(1 -1,1 4)") as LineString;
            var polygon = reader.Read("Polygon((0 0,0 1, 1 1,1 0,0 0))") as Polygon;
        }
    }
}

注意:

  • Read里面字符串中的"Point"、“LineString”、"Polygon"不区分大小写。

读取shp文件

将shp数据读取到FeatureCollection(可以转到定义之处查看常用方法)对象中,再进行后续操作。
NTS1.14.0以上的版本似乎没有FeatureCollection(1.15.0版本中发现没有)。

public static FeatureCollection ReadShpFile(string pathName)
{
    FeatureCollection featureCollection = new FeatureCollection();
    IGeometryFactory gfactory = GeometryFactory.Floating;
    ShapefileDataReader dataReader = new ShapefileDataReader(pathName, gfactory);

    while (dataReader.Read())
    {
        Feature feature = new Feature { Geometry = dataReader.Geometry };
        int length = dataReader.DbaseHeader.NumFields;
        var keys = new string[length];
        for (var i = 0; i < length; i++)
        {
            keys[i] = dataReader.DbaseHeader.Fields[i].Name;
        }
        feature.Attributes = new AttributesTable();
        for (var i = 0; i < length; i++)
        {
            var val = dataReader.GetValue(i);
            string value = val.ToString();
            feature.Attributes.AddAttribute(keys[i], value);
        }
        featureCollection.Add(feature);
    }
    return featureCollection;
}

生成shp文件

private static void OutputShp(string shpPath, List<IGeometry> PolygonList)
{
    //要素集
    FeatureCollection featureCollection = new FeatureCollection();

    //遍历
    string wkt;
    foreach (var polygon in PolygonList)
    {
        //根据wkt制作一个几何
        wkt = polygon.AsText();
        var geometry = wktReader.Read(wkt);

        //创建一个要素
        var feature = new Feature();
        featureCollection.Add(feature);

        //设置几何图形
        feature.Geometry = geometry;

        //设置属性表
        var attrTable = new AttributesTable();
        feature.Attributes = attrTable;//属性表哪怕是空的,也不能缺少
        //默认有两个属性列,至少再多添加一个属性,比如id。
        //这样才能较好地在ArcMap中使用该数据时(比如编辑时不会卡住)
        attrTable.AddAttribute("id", "");
    }

    IGeometryFactory geometryFactory = GeometryFactory.Floating;//Floating起什么作用?
    var dataWriter = new ShapefileDataWriter(shpPath, geometryFactory, Encoding.Default);//设置Encoding编码方式
    dataWriter.Header = ShapefileDataWriter.GetHeader(featureCollection[0], featureCollection.Count);
    dataWriter.Write(featureCollection.Features);
}

几何相关

线段上内插点

点集的顺时针或逆时针的判断

public Coordinate Interpolation(LineSegment segment, double dist)
{
	//按dist距离,从起点开始内插得到一个点,返回这个点
    return segment.PointAlong(dist / segment.Length);
}
Coordinate[] coords = new Coordinate[]
{
    new Coordinate(0,0),
    new Coordinate(10,0),
    new Coordinate(20,0),
    new Coordinate(10,0),
    new Coordinate(10,10),
    new Coordinate(0,10),
    new Coordinate(0,0),//首尾必须完成闭合,才能成功创建线性环
};
var ring = new LinearRing(coords);//线性环
var isCounterClockwise = ring.IsCCW;//true,逆时针。CCW是CounterClockwise的缩写,表示逆时针

矢量相关

转角方式介绍

  • 右手准则
  • 逆时针为正
  • 顺时针为负
  • 单位弧度
Vector2D vec1 = new Vector2D(1, 0);
Vector2D vec2 = new Vector2D(0, -1);

var angle12 = vec1.AngleTo(vec2);
var angle21 = vec2.AngleTo(vec1);

在这里插入图片描述

其他

var v2 = vec1 .RotateByQuarterCircle(1);//逆时针旋转一个四分之一圆,即90度
var v3 = vec1 .Rotate(Math.PI / 2);//逆时针旋转90度,必须是弧度
var v3 = v2.Normalize();//获得单位矢量
var v4 = v3.Multiply(2.5);//乘上倍数
var v = new Vector2D(1,0);
var v1 = new Vector2D(1,1);
var v2 = new Vector2D(0,1);
var v3 = new Vector2D(-1,1);
var v4 = new Vector2D(-1,0);
var v5 = new Vector2D(-1,-1);
var v6 = new Vector2D(0, -1);
var v7 = new Vector2D(1,-1);

var c1 = v.AngleTo(v1);//π/4
var c2 = v.AngleTo(v2);//π/2
var c3 = v.AngleTo(v3);//3π/4
var c4 = v.AngleTo(v4);//π
var c5 = v.AngleTo(v5);//-3π/4
var c6 = v.AngleTo(v6);//-π/2
var c7 = v.AngleTo(v7);//-π/4

算法相关

泰森多边形

介绍:

  • 泰森多边形,即Voronoi图

  • NTS提供了一个工具类VoronoiDiagramBuilderAPI链接

  • 输入一些点,就可以创建泰森多边形。

  • 返回的Polygon的集合

例子:

//Voronoi图
VoronoiDiagramBuilder vor = new VoronoiDiagramBuilder();

//生成的范围(如果太小,将不起作用,而是默认大小)
var enve = new Envelope(new Coordinate(0, 0), new Coordinate(15, 15));
vor.ClipEnvelope = enve;

//添加点数据
List<Coordinate> pts = new List<Coordinate>()
{
    new Coordinate(0,0),
    new Coordinate(10,0),
    new Coordinate(5,10),
};
vor.SetSites(pts);

//生成结果           
var mp = vor.GetDiagram(GeometryFactory.Floating);

//遍历结果,根据需要处理每一个多边形(Polygon)
List<Polygon> geos = new List<Polygon>();
foreach (var geo in mp.Geometries)
{
    var polygon = geo as Polygon;
    geos.Add(polygon);

    //可以获得该多边形中的点
    //(也就是一开始通过SetSite方法,参与创建多边形的点,获得的点与原数据只是坐标数值相同,不具有相同的引用)
    var coord = geo.UserData as Coordinate;
}

利用前面的生成shp的方法,将得到的List<Polygon> geos ,也就是多边形数据输出shp文件,并在ArcMap软件中浏览(绿色点是额外添加的):
在这里插入图片描述
应用:参考博客实现C#版的多边形中心线的提取。

道格拉斯普克简化面

道格拉斯普克法简化几何图形。

比如简化一个多边形:

var polygonNew = DouglasPeuckerSimplifier.Simplify(polygon, 20);//20为距离容差

空间运算相关

线与线相交、相减、对称差、合并

1. 十字型交点
在这里插入图片描述

WKTReader reader = new WKTReader();
var L1 = reader.Read("LineString(0 1,3 1)") as LineString;
var L2 = reader.Read("LineString(1 0,1 3)") as LineString;
var r1 = L1.Intersection(L2);//相交
var r2 = L1.Difference(L2);//相减
var r3 = L1.SymmetricDifference(L2);//对称差
var r4 = L1.Union(L2);//合并

断点处看参数:
在这里插入图片描述

  • Intersection:相交,返回一个点坐标。(交集)
  • Difference:相减,有先后顺序之别。上面是L1减去L2,即L1减去其与L2的交点,得到两条线段。(单方减去交集)
  • SymmetricDifference:对称差,L1+L2,再减去两者的交点,得到4条线段。(并集减去交集)
  • Union:合并,L1+L2,因为有交点,所以分成了4条线段。(并集)

2. L型交点
在这里插入图片描述
其余不变,就改一下点坐标:

var L1 = reader.Read("LineString(1 1,3 1)") as LineString;
var L2 = reader.Read("LineString(1 1,1 3)") as LineString;

L型和十字型的结果,所反映的情况分析是一样的:

在这里插入图片描述
可以看到Union这个方法并不能将两条线合并为一条,即使它们有共同的邻点。

3. 无交点
在这里插入图片描述
继续修改点坐标:

var L1 = reader.Read("LineString(1 1,3 1)") as LineString;
var L2 = reader.Read("LineString(1 2,1 3)") as LineString;

由于没有交点,所以Intersection结果为EMPTY(空)。
相减、对称差、合并则在此空的基础上进行。

在这里插入图片描述

面与面相交、相减、对称差、合并

  1. 相交
    如下图,面S1和面S2相交,交集为S3:
    在这里插入图片描述

    WKTReader reader = new WKTReader();//用于读取wkt格式的几何字符串
    var S1 = reader.Read("Polygon((0 0,0 2,4 2,4 0,0 0))") as Polygon;//将字符串解析为几何对象
    var S2 = reader.Read("Polygon((1 1,1 3,2 3,2 1,1 1))") as Polygon;
    var r1 = S1.Intersection(S2);//相交
    var r2 = S1.Difference(S2);//相减:S1-S2,前面的减去后面的
    var r3 = S1.SymmetricDifference(S2);
    var r4 = S1.Union(S2);
    

    在这里插入图片描述

    • Intersection:相交,得到S3,就是S1与S2的交集。
    • Difference:相减,S1-S2=S1-S3,就是S1去掉S3变凹的图形。
    • Symmetric:对称差,交集取反。S1+S2得到一个凸的图形,再减去S3,变成两个图形。
    • Union:合并,交集,得到一个凸的图形。

2. 相切

修改点坐标:

var S1 = reader.Read("Polygon((0 0,0 2,4 2,4 0,0 0))") as Polygon;
var S2 = reader.Read("Polygon((1 2,1 4,2 4,2 2,1 2))") as Polygon;

S1和S2相切,有一条共边L1。
在这里插入图片描述
在这里插入图片描述

  • 相交,得到一条交线。
  • 相减,S1-S2=S1-L1=S1,减去交线并不影响S1,结果还是S1。
  • 对称差,S1+S2得到合并后的凸的图形,再减去交线L1,也就没有变化。
  • 合并,S1+S2,和对称差的结果相同。

3. 相离

修改点坐标:

var S1 = reader.Read("Polygon((0 0,0 2,4 2,4 0,0 0))") as Polygon;
var S2 = reader.Read("Polygon((1 3,1 5,2 5,2 3,1 3))") as Polygon;

S1和S2相离:
在这里插入图片描述
在这里插入图片描述

  • 相交为空。
  • 相减,S1-S2=S1-空=S1,就是S1本身。
  • 对称差,S1+S2-交集=S1+S2-空集=S1+S2,得到两个面,只是简单地放到一个集合中(也就是属于MULTIPOLYGON类型的面),本身无法合并。
  • 合并,S1+S2,得到两个面(也是MULTIPOLYGON类型)。

二维几何比较

Coordinate

Coordinate c1 = new Coordinate(10.01, 1.01);
Coordinate c2 = new Coordinate(10.01, 1.01);

var aa = c1.Equals2D(c2);//true
var bb = c1.Equals(c2);//true
var cc = c1 == c2;//false

LineSegment

// 起点和终点的顺序可以不同
private bool CoordEquals(LineSegment L1, LineSegment L2)
{
    if (L1.P0.Equals2D(L2.P0) && L1.P1.Equals2D(L2.P1))
    {
        return true;
    }
    if (L1.P0.Equals2D(L2.P1) && L1.P1.Equals2D(L2.P0))
    {
        return true;
    }
    return false;
}
private bool CoordEquals(LineSegment L1, Coordinate c1, Coordinate c2)
{
    if (L1.P0.Equals2D(c1) && L1.P1.Equals2D(c2))
    {
        return true;
    }
    if (L1.P0.Equals2D(c2) && L1.P1.Equals2D(c1))
    {
        return true;
    }
    return false;
}

计算速度粗略比较

Coordinate和Point的Distance的比较

Coordinate c1 = new Coordinate(10, 20);
Coordinate c2 = new Coordinate(103, 2333);
Point p1 = new Point(c1);
Point p2 = new Point(c2);

MyTimer.Prepare();
MyTimer.Start();
for (int i = 0; i < 100000; i++)
{
	var a = c1.Distance(c2);
}
MyTimer.PrintDuration();
MyTimer.Start();
for (int i = 0; i < 100000; i++)
{
	var a = p1.Distance(p2);
}
MyTimer.PrintDuration();

打印的结果是:

0.000356
0.0465025

就上述情况而言,用Coordinate对象计算距离更合适。

Polygon包含LineString、Point的速度比较

Coordinate[] coords = new Coordinate[]
{
	new Coordinate(0,0),
	new Coordinate(0,10),
	new Coordinate(10,10),
	new Coordinate(10,7),
	new Coordinate(3,7),
	new Coordinate(3,0),
	new Coordinate(0,0),
};
var ring = new LinearRing(coords);
var polygon2 = new Polygon(ring);

Coordinate c1 = new Coordinate(2, 9);
Coordinate c2 = new Coordinate(6, 8);
Coordinate c3 = new Coordinate(16, 8);
Point p1 = new Point(c1);
Point p2 = new Point(c2);
Point p3 = new Point(c3);

LineString line1 = new LineString(new Coordinate[] { c1, c2 });
LineString line2 = new LineString(new Coordinate[] { c1, c3 });

MyTimer.Prepare();
MyTimer.Start();
for (int i = 0; i < 100000; i++)
{
	polygon2.Contains(p1);
	polygon2.Contains(p2);
	polygon2.Contains(p1);
	polygon2.Contains(p3);
}
MyTimer.PrintDuration();
MyTimer.Start();
for (int i = 0; i < 100000; i++)
{
	polygon2.Contains(line1);
	polygon2.Contains(line2);
}
MyTimer.PrintDuration();

打印结果:

4.9734832
2.5062652

试验的数据量、数据分布情况较简单。就上述情况而言,用LineString被包含在面中反而比它的两个端点Point被包含在面中计算得更快。

缓冲区

默认的缓冲区参数(圆角)

var S = GeoMgr.Read("polygon((0 0,0 5,1 5,1 1,4 1,4 5,5 5,5 0,0 0))") as Polygon;
var T = GeoMgr.Read("polygon((0 -5,0 5,2 5, 3 0, 2 -5,0 -5))") as Polygon;
var L = GeoMgr.Read("linestring(1 5, 4 5)") as LineString;

var S2 = S.Buffer(0.1);
var L2 = L.Buffer(0.1);
var T2 = T.Buffer(0.1);

效果图:

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

自定义缓冲区参数(斜角)

//数据
var S = GeoMgr.Read("polygon((0 0,0 5,1 5,1 1,4 1,4 5,5 5,5 0,0 0))") as Polygon;
var T = GeoMgr.Read("polygon((0 -5,0 5,2 5, 3 0, 2 -5,0 -5))") as Polygon;
var L = GeoMgr.Read("linestring(1 5, 4 5)") as LineString;

//缓冲区参数
var bufferParam = new BufferParameters();
bufferParam.JoinStyle = JoinStyle.Mitre;
bufferParam.EndCapStyle = EndCapStyle.Flat;

//缓冲区结果
var S2 = S.Buffer(-0.1, bufferParam);
var L2 = L.Buffer(-0.1, bufferParam);
var T2 = T.Buffer(-0.1, bufferParam);

效果图:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

补充

缓冲区距离可以为负值,这样会向内做缓冲区。

评论 29
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值