1 using System.Threading;
2 using System.Collections.Generic;
3 using GeoAPI.Geometries;
4 using NetTopologySuite.Geometries;
5 using NetTopologySuite.Operation.Polygonize;
6 using System.Linq;
7
8 namespace Test.Workspace
9 {
10 public static class NTSExpend
11 {
12 public static IGeometry Validate(this IGeometry geometry) => NTSUtils.Instance.Validate(geometry);
13 }
14 public class NTSUtils
15 {
16 private static NTSUtils _this = null;
17 private static object lockObj = new object();
18
19 /// <summary>
20 /// 单例创建对象
21 /// </summary>
22 /// <returns></returns>
23 public static NTSUtils Instance
24 {
25 get
26 {
27 if (_this == null)
28 {
29 lock (lockObj)
30 {
31 if (_this == null)
32 _this = new NTSUtils();
33 return _this;
34 }
35 }
36 else
37 return _this;
38 }
39 }
40
41 /// <summary>
42 /// 处理图形自相交
43 /// </summary>
44 /// <returns></returns>
45 public IGeometry Validate(IGeometry geometry)
46 {
47 Polygonizer polygonizer = null;
48 if (geometry is Polygon)
49 {
50 if (geometry.IsValid)
51 {
52 geometry.Normalize();
53 return geometry;
54 }
55 polygonizer = new Polygonizer();
56
57 AddPolygon((Polygon)geometry, polygonizer);
58 }
59 else if (geometry is MultiPolygon)
60 {
61 if (geometry.IsValid)
62 {
63 geometry.Normalize();
64 return geometry;
65 }
66 polygonizer = new Polygonizer();
67 for (int n = geometry.NumGeometries; n-- > 0;)
68 {
69 AddPolygon((Polygon)geometry.GetGeometryN(n), polygonizer);
70 }
71 }
72 else if (geometry is GeometryCollection)
73 {
74 if (geometry.IsValid)
75 {
76 geometry.Normalize();
77 return geometry;
78 }
79 polygonizer = new Polygonizer();
80 for (int n = geometry.NumGeometries; n-- > 0;)
81 {
82 IGeometry temp = geometry.GetGeometryN(n);
83 if (temp is Polygon)
84 {
85 AddPolygon((Polygon)temp, polygonizer);
86 }
87 else if (temp is MultiPolygon)
88 {
89 for (int m = temp.NumGeometries; m-- > 0;)
90 {
91 AddPolygon((Polygon)temp.GetGeometryN(m), polygonizer);
92 }
93 }
94 }
95 }
96 else
97 return geometry;
98
99 return ToPolygonGeometry(polygonizer.GetPolygons(), geometry.Factory);
100 }
101
102 void AddPolygon(Polygon polygon, Polygonizer polygonizer)
103 {
104 AddLineString(polygon.ExteriorRing, polygonizer);
105 for (int n = polygon.NumInteriorRings; n-- > 0;)
106 {
107 AddLineString(polygon.GetInteriorRingN(n), polygonizer);
108 }
109 }
110
111 void AddLineString(ILineString lineString, Polygonizer polygonizer)
112 {
113
114 if (lineString is LinearRing)
115 { // LinearRings are treated differently to line strings : we need a LineString NOT a LinearRing
116 lineString = lineString.Factory.CreateLineString(lineString.CoordinateSequence);
117 }
118
119 // unioning the linestring with the point makes any self intersections explicit.
120
121 IPoint point = lineString.Factory.CreatePoint(lineString.GetCoordinateN(0));
122 IGeometry toAdd = lineString.Union(point);
123
124 //Add result to polygonizer
125 polygonizer.Add(toAdd);
126 }
127
128 /**
129 * Get a geometry from a collection of polygons.
130 *
131 * @param polygons collection
132 * @param factory factory to generate MultiPolygon if required
133 * @return null if there were no polygons, the polygon if there was only one, or a MultiPolygon containing all polygons otherwise
134 */
135 IGeometry ToPolygonGeometry(ICollection<IGeometry> polygons, IGeometryFactory factory)// geoSource.factory 可以创建多边形
136 {
137 switch (polygons.Count())
138 {
139 case 0:
140 return null; // No valid polygons!
141 case 1:
142 return polygons.FirstOrDefault(); // single polygon - no need to wrap
143 default:
144 //polygons may still overlap! Need to sym difference them
145 List<IGeometry> ps = polygons.ToList();
146 double maxArea = ps.Select(x => x.Envelope.Area).Max();
147 IGeometry geometry = ps.FirstOrDefault(x => x.Envelope.Area == maxArea);
148 IGeometry result = geometry;
149 for (int i = 0; i < ps.Count; i++)
150 {
151 IGeometry currentGeo = ps[i];
152
153 //string testStr = currentGeo.ToString();
154 //string res = GdalShapeOperate.Instance.GetPoints(OSGeo.OGR.Ogr.CreateGeometryFromWkt(ref testStr, null));
155 if (currentGeo == null || currentGeo.Equals(geometry))
156 continue;
157
158 //1. 判断当前图形是否是多环图形
159
160 //Covers 优于 Contains
161 // symDifference 并集去掉重合部分
162 IPolygon p = result as IPolygon;
163 if (p?.NumInteriorRings > 0)
164 {
165 //多环的话 需要把内部空环依次和当前图形比对 判断是否是重叠 (重叠取diff)
166 Polygonizer polygonizer = new Polygonizer();
167 for (int n = p.NumInteriorRings; n-- > 0;)
168 {
169 AddLineString(p.GetInteriorRingN(n), polygonizer);
170 }
171 List<IGeometry> interPolygons = polygonizer.GetPolygons().ToList();
172 //内环存在的话 需要和当前图形比对 如果基本匹配(重合面积>90) 取差集
173 if (interPolygons.Count > 0)
174 {
175 bool isIntersect = false;
176 foreach (var item in interPolygons)
177 {
178 if (item.Intersects(currentGeo) && item.Intersection(currentGeo).Area / currentGeo.Area > 0.9)
179 {
180 isIntersect = true;
181 break;
182 }
183 }
184 if (isIntersect)
185 {
186 result = result.Difference(currentGeo);
187 continue;
188 }
189 }
190 }
191
192 //不是多环图形, 如果重叠 且重叠面积>当前图形面积的90%则 交集, 否则并集
193 if (geometry.Intersects(currentGeo) && geometry.Intersection(currentGeo).Area / currentGeo.Area > 0.9)
194 result = result.Difference(currentGeo);
195 else
196 result = result.Union(currentGeo);
197
198 }
199 return result;
200 }
201 }
202
203 /// <summary>
204 /// wkt转geometry
205 /// </summary>
206 /// <param name="wkt"></param>
207 /// <returns></returns>
208 public IGeometry WktToGeometry(string wkt)
209 {
210 NetTopologySuite.IO.WKTReader reader = new NetTopologySuite.IO.WKTReader();
211 IGeometry geo = reader.Read(wkt);
212 return geo;
213 }
214
215 /// <summary>
216 /// rings 转 geometry
217 /// </summary>
218 /// <param name="rings"></param>
219 /// <returns></returns>
220 public IGeometry RingsToGeometry(string rings)
221 {
222 OSGeo.OGR.Geometry geometry = GdalShapeOperate.Instance.GetGeometryFromShapeJson(rings);
223
224 string wkt = string.Empty;
225 while (wkt == string.Empty)
226 {
227 wkt = GdalShapeOperate.Instance.ExportToWkt(geometry);
228 Thread.Sleep(50);
229 }
230 return WktToGeometry(wkt);
231 }
232
233 public string WktToRings(string wkt) => GdalShapeOperate.Instance.GetPoints(OSGeo.OGR.Geometry.CreateFromWkt(wkt));
234
235 public string GeometryToRings(Geometry geo) => GdalShapeOperate.Instance.GetPoints(OSGeo.OGR.Geometry.CreateFromWkt(geo.ToString()));
236 }
237 }
