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:

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;

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:

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;

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.

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);

        }
    }
}

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.

Published by webfresh