Spikes: research with ESRI

A short follow-up on my previous blog about spikes and their causes. My collegue Bert recoded the scenario in C# using ESRI. Thanks Bert! The full code is shown below.

The answer is this: with ESRI the spikes are not there. The program produces:

area: 295.44447138901, perimeter: 92.5058167188044

select geometry::STGeomFromText ('POLYGON((180956.436999999
313716.998,180953.978700001 313701.0152,180926.420299999
313732.510699999,180954.359000001
313720.504000001,180956.436999999 313716.998,180956.436999999
313716.998))', 28992)

And this result is right, area and perimeter are without spikes and if you visualize this WKT (e.g. with SQL Server) you get the right picture...

Interesting is that all coordinates seems to be rounded on one millimeter or a tenth of it. That is probably the result of the fuzzy tolerance policy that ESRI uses. Personally I'm normally a bit sceptic about that approach: though it can avoid spikes and slivers, it also brings some fuzzyness into the result. It is sometimes quite hard to determine a tolerance which is sufficient for cleaning and sufficient for not destroying features... Anyway, in this scenario it works well.

It seems interesting to check, if you use rounding in e.g. Boost.Geometry, would the spikes be gone... The answer is (of course), sometimes yes, sometimes no. The rounding is applied on a larger grid, so the effect remains the same. Rounding is not enough. To get rid of spikes, in this scenario, the spike removal functionality is advised. Which uses a gap distance, which can be seen as a sort of fuzzy tolerance...

C# Code for "difference after cut" with ESRI

using System; 
using System.Collections.Generic; 
using System.ComponentModel; 
using System.Data; 
using System.Drawing; 
using System.Text; 
using System.Windows.Forms; 
using ESRI.ArcGIS.Geometry; 
using ESRI.ArcGIS.esriSystem; 
using Microsoft.SqlServer.Types; 
using System.Data.SqlTypes; 

namespace SpikesResearchWithEsri 
{ 
    public partial class Form1 : Form
    { 
        public Form1() 
        { 
          InitializeComponent(); 
        } 

        private void InitializeLicense()
        { 
          AoInitialize aoi = new AoInitializeClass(); 

          //Additional license choices can be included here. 
          esriLicenseProductCode productCode = 
          esriLicenseProductCode.esriLicenseProductCodeArcEditor; 
          if (aoi.IsProductCodeAvailable(productCode) == 
          esriLicenseStatus.esriLicenseAvailable) 
          { 
          aoi.Initialize(productCode); 
          } 
        } 


        private void button1_Click(object sender, EventArgs e)
        { 


          ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop); 
          InitializeLicense(); 
          IPolygon parcel = this.GetParcel(); 
          IPolyline line = this.GetCutLine(); 

          IPolygon bufferedLine=(IPolygon)((ITopologicalOperator)line).Buffer(0.1); 
          
          IGeometry result=((ITopologicalOperator)parcel).Difference(bufferedLine); 
          IPolygon4 p = (IPolygon4)result; 

          // Big part: 
          IPolygon firstPolygon = this.GetPolygon(p, true);
          // Small part: 
          IPolygon secondPolygon = this.GetPolygon(p, false);


          // use if we do a cutpolygon 
          //IGeometry leftPolygon; 
          //IGeometry rightPolygon; 
          //((ITopologicalOperator)parcel).Cut(line, out leftPolygon, out rightPolygon); 
          //IPolyline cutline=(IPolyline)((ITopologicalOperator)(ICurve)parcel).Intersect(line, esriGeometryDimension.esriGeometryNoDimension); 
          //IPolygon firstPolygon = (IPolygon)leftPolygon; 
          //IPolygon secondPolygon = (IPolygon)rightPolygon; 
          //this.textBox1.Text=("from: " + cutline.FromPoint.X.ToString() + ", " + 
          // cutline.FromPoint.Y.ToString() + ", " + 
          // "to: " + cutline.ToPoint.X.ToString() + ", " + 
          // cutline.ToPoint.Y.ToString()) + Environment.NewLine; 

          this.textBox1.Text += "Small polygon stats: " + System.Environment.NewLine; 
          this.textBox1.Text += GetStats(secondPolygon); 
          this.textBox1.Text += GetWkt(secondPolygon); 

          // now calculate the difference between original and the big part (result should be the small part) 
          IPolygon4 pol2= (IPolygon4)((ITopologicalOperator)parcel).Difference(firstPolygon); 
          this.textBox1.Text += "Result Small polygon stats: " + System.Environment.NewLine; 
          this.textBox1.Text += GetStats(pol2); 
          this.textBox1.Text += GetWkt(pol2); 
        } 

        private string GetStats(IPolygon pol)
        { 
          double area = ((IArea)pol).Area; 
          double perimeter = ((ICurve)pol).Length; 
          return "area: " + area.ToString() + ", perimeter: " + perimeter.ToString() + System.Environment.NewLine; 
        } 


        private string GetWkt(IPolygon pol)
        { 
          string res = GeomHelper.aoGeomToWkt(pol, false, 28992);
          res = res.Substring(0, res.Length - 1);
          string print = "select geometry::STGeomFromText ('" + res + "', 28992)" + Environment.NewLine; 
          return print; 
        } 

        public IPolygon GetPolygon(IPolygon4 input, bool first)
        { 
          IGeometryBag exteriorRings = input.ExteriorRingBag; 
          IEnumGeometry exteriorRingsEnum = exteriorRings as IEnumGeometry;
          exteriorRingsEnum.Reset(); 

          IRing firstRing = exteriorRingsEnum.Next() as IRing;
          IPointCollection pcFirst = (IPointCollection)firstRing; 

          IRing secondRing = exteriorRingsEnum.Next() as IRing;
          IPointCollection pcSecond = (IPointCollection)secondRing; 

          IPolygon p = new PolygonClass(); 

          if(first)((IPointCollection)p).AddPointCollection(pcFirst); 
          if(!first)((IPointCollection)p).AddPointCollection(pcSecond); 
          return p; 
        } 


        private IPolygon GetParcel() 
        { 
          string wkt = "POLYGON ((180956.437 313716.998, 180954.359 313720.504, 180906.303 313741.156, 180896.82 313744.272, 180890.456 313748.17, 180831.49 313797.135, 180809.668 313816.098, 180776.418 313842.336, 180771.483 313846.102, 180746.934 313843.764, 180731.738 313840.386, 180717.32 313833.763, 180710.827 313828.438, 180704.982 313820.515, 180699.267 313811.813, 180698.488 313801.162, 180700.826 313793.889, 180705.539 313791.06, 180706.021 313790.771, 180710.156 313789.208, 180766.157 313768.042, 180825.515 313744.013, 180878.636 313724.271, 180945.593 313695.884, 180952.801 313693.358, 180956.437 313716.998))"; 
          IGeometryFactory factory = new GeometryEnvironmentClass();
          SqlGeometry g = SqlGeometry.Parse(wkt); 
          IGeometry geom2 = new PolygonClass(); 
          int countout; 
          factory.CreateGeometryFromWkbVariant(g.STAsBinary().Value, out geom2, out countout);
          geom2.SpatialReference = this.GetProjectedSpatialReference(28992);
          return (IPolygon) geom2; 
        } 

        private IPolyline GetCutLine() 
        { 
          string line = "LINESTRING(180955 313700,180920 313740)"; 
          IGeometryFactory factory = new GeometryEnvironmentClass();
          SqlGeometry g = SqlGeometry.Parse(line); 
          IGeometry geom2 = new PolylineClass(); 
          geom2.SpatialReference = this.GetProjectedSpatialReference(28992);
          int countout; 
          
          factory.CreateGeometryFromWkbVariant(g.STAsBinary().Value, out geom2, out countout);
          return (IPolyline)geom2; 

        } 


        private ISpatialReference GetProjectedSpatialReference(int pcsType)
        { 
          ISpatialReferenceFactory pSRF = new SpatialReferenceEnvironmentClass();
          IProjectedCoordinateSystem m_ProjectedCoordinateSystem = pSRF.CreateProjectedCoordinateSystem(pcsType);
          ISpatialReference spatialReference = (ISpatialReference)m_ProjectedCoordinateSystem; 
          return spatialReference; 
        } 

    } 
}