文章目录
介绍
简称NTS，是C#端的GIS相关的拓扑库。Github地址、API文档。
其他相关
关于OGC标准（开放地理空间信息联盟(Open Geospatial Consortium)）
项目环境
- 本文以VS2019创建的.Net Framework4.7.2 控制台应用为例，对NTS的一些功能使用进行介绍。
- 右击项目的引用，打开NuGet程序包管理器，搜索
NetTopologySuite，安装NetTopologySuite1.14.0版和NetTopologySuite.IO1.14.0.1版本。
如果要看源码，考虑到网速问题，可以去Gitee上（已从Github中迁移过来，访问地址），根据分支中的标签，选择与NuGet中相同版本的源码。
常用的API介绍
IO相关
读取wkt字符串
作用：WKTReader将WKT格式的字符串转化为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图
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（空）。
相减、对称差、合并则在此空的基础上进行。
面与面相交、相减、对称差、合并
-
相交
如下图，面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);
效果图：
补充
缓冲区距离可以为负值，这样会向内做缓冲区。
扫一扫
939
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
