Simplifying Spatial Data while Preserving Topologies in C#
The data in a shape file is frequently too detailed for use over a web page / internet connection, as the amount of data being transferred is large, but the user typically views the data on a zoomed-out map, and isn’t in need of all the detail. MapShaper is a fantastic tool for converting shape files to other formats, while simplifying the shapes in the process. However, it’s either a manual process using their website, or there is a command line tool written in JavaScript. But can you do it in C#/.net?
One of the main challenges is to simplify the shapes together. Each shape in the shape file is a different record/row, and if you simplify them individually, the shared boundaries between the shapes will no longer be shared, and you will end up with overlapping boundaries or gaps.
Take the following code example:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
var entities = await db.Entities.ToListAsync(); dynamic geoJson = new ExpandoObject(); geoJson.type = "FeatureCollection"; var features = new List<dynamic>(); foreach (var entity in entities) { // simplify (uses SqlGeometry.Reduce) var simplified = Utilities.Geo.Simplify(entity.Boundaries, tolerance); // convert to GeoJSON4EntityFramework.Feature with the appropriate Id var feature = Utilities.Geo.GetFeature(simplified, entity.Id); features.Add(feature); } geoJson.features = features; return geoJson; |
The resulting GeoJson looks like the following:
If you Union (join) the shapes before simplifying, such as in this code:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
var entities = await db.Entities.ToListAsync(); // bit of a hack to union all the boundaries DbGeometry allBoundaries = null; for (var i = 0; i < entities.Count; i++) { if (i == 0) allBoundaries = entities[i].Boundaries; else allBoundaries = allBoundaries.Union(entities[i].Boundaries); } // simplify (uses SqlGeometry.Reduce) var simplified = Utilities.Geo.Simplify(allBoundaries, tolerance); dynamic geoJson = new ExpandoObject(); geoJson.type = "FeatureCollection"; var features = new List<dynamic>(); // convert to GeoJSON4EntityFramework.Feature with the (in)appropriate Id var feature = Utilities.Geo.GetFeature(simplified, "ALL"); features.Add(feature); geoJson.features = features; return geoJson; |
Your GeoJson will be produced as follows:
The Union function actually merges the shapes, so you end up with a single overall shape that is then simplified.
So, how do you get the required result?
One approach is to convert the polygons and multipolygons into lines (LineStrings and MultiLineStrings), then simplify the lines (while preserving the topologies). Then you have to convert the areas between the lines back into polygons, and then re-attach any properties (e.g. the Id or Name of the shape).
The first few codes samples above made use of the spatial functions directly in SQL Server. However, a limitation in SQL Server is that you can’t (easily) convert from the simplified lines back to polygons. So our approach was to use NetTopologySuite, a C# library that is a direct port from Java Topology Suite, which offers a bit more functionality than what you get natively in SQL Server (using SqlGeometry etc.).
Here is a simplified(!) version of the code, in a console app:
Note that there are 3 items you’d need to replace:
- YOUR_USERNAME: your windows username, or change the path to something suitable.
- DC31_COORDINATES & DC32_COORDINATES: this demo just uses two shapes/polygons: you can paste your co-ordinates into these placeholders, so that it uses inline well-known-text values, or you can populate the polygons in another way, e.g. from a database field, or a shape file.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 |
using GeoAPI.Geometries; using NetTopologySuite.Features; using NetTopologySuite.Geometries; using NetTopologySuite.IO; using NetTopologySuite.Operation.Linemerge; using NetTopologySuite.Operation.Polygonize; using NetTopologySuite.Precision; using NetTopologySuite.Simplify; using System; using System.Collections.Generic; using System.IO; using System.Linq; namespace Simplifier { class Program { public const string OutputFolder = @"C:\Users\YOUR_USERNAME\Desktop\"; static void Main(string[] args) { // will hold all the polygons var polygons = new Dictionary<string, Polygon>(); var factory = new OgcCompliantGeometryFactory(); var wktReader = new WKTReader(factory); // load from wkt polygons.Add("DC31", (Polygon)wktReader.Read("POLYGON((DC31_COORDINATES))")); polygons.Add("DC32", (Polygon)wktReader.Read("POLYGON((DC32_COORDINATES))")); // get exterior rings and merge them into a single geometry var unionedExteriorRings = polygons["DC31"].ExteriorRing.Union(polygons["DC32"].ExteriorRing); OutputAsGeoJson(unionedExteriorRings, "unionedExteriorRings"); // merge the line strings var merger = new LineMerger(); merger.Add(unionedExteriorRings); var mergedLineStrings = merger.GetMergedLineStrings(); OutputAsGeoJson(mergedLineStrings, "mergedLineStrings"); var geometryCollection = new GeometryCollection(mergedLineStrings.ToArray()); // simplify the linestrings (exterior rings) var simplified = Simplify(geometryCollection); OutputAsGeoJson(simplified, "simplified"); // convert back to polygons var polygonizer = new Polygonizer(false); polygonizer.Add(simplified); var newPolygons = polygonizer.GetPolygons(); OutputAsGeoJson(newPolygons, "newPolygons"); // make valid if not right-hand-rule // nb: currently only works for polygons var valid = newPolygons.Select(o => ((Polygon)o).Shell.IsCCW ? o : o.Reverse()); OutputAsGeoJson(valid, "valid"); // reduce precision as i don't need 15 digits var reduced = valid.Select(o => Reduce(o)); OutputAsGeoJson(reduced, "reduced"); // match back to the original table using the area of intersection // if it intersects by more than 50%, then assume it's the same entity var results = new List<GeometryWithId>(); foreach (var r in reduced) { var result = new GeometryWithId { Geometry = r }; foreach (var poly in polygons) { if (poly.Value.Intersects(result.Geometry) && poly.Value.Intersection(result.Geometry).Area / result.Geometry.Area > 0.5 ) result.Id = poly.Key; } results.Add(result); } OutputAsGeoJson(results, "results"); Console.WriteLine("finished...."); Console.ReadKey(); } private class GeometryWithId { public string Id; public IGeometry Geometry; } static void OutputAsGeoJson(IEnumerable<GeometryWithId> results, string fileName) { var geoJson = new GeoJsonWriter(); var features = new FeatureCollection(); foreach (var result in results) { var feature = new Feature(); feature.Geometry = result.Geometry; var attributesTable = new AttributesTable(); attributesTable.Add("Id", result.Id); feature.Attributes = attributesTable; features.Add(feature); } var json = geoJson.Write(features); File.WriteAllText(OutputFolder + $@"{fileName}.json", json); } static void OutputAsGeoJson(IGeometry geometry, string fileName) { OutputAsGeoJson(new List<IGeometry> { geometry }, fileName); } static void OutputAsGeoJson(IEnumerable<IGeometry> geometries, string fileName) { OutputAsGeoJson(geometries.Select(o => new GeometryWithId { Geometry = o }), fileName); } static IGeometry Simplify(IGeometry geometry) { /// Note that in general D-P does not preserve topology - /// e.g. polygons can be split, collapse to lines or disappear /// holes can be created or disappear, /// and lines can cross. /// To simplify point while preserving topology use TopologySafeSimplifier. /// (However, using D-P is significantly faster). //return NetTopologySuite.Simplify.DouglasPeuckerSimplifier.Simplify(geometry, 0.1); return TopologyPreservingSimplifier.Simplify(geometry, 0.1); } static IGeometry Reduce(IGeometry geometry) { // keep only 4 decimal places return new GeometryPrecisionReducer(new PrecisionModel(1e4)).Reduce(geometry); } } } |
The code will generate GeoJson files at multiple points in the processing. The last step produces results.json, which can be viewed in geojsonlint.com, as per the screenshot below.
You can see the two shapes (in this case districts DC31 & DC32) have been simplified, while the shared edge has been preserved (no overlaps or gaps), and the attributes (id) have also been retained.
The above code will need a bit of tweaking for more complex geometries, but it demonstrates the basic premise.
This webpage provided the basic solution in PostGIS.