Buffer playground

buffer_playground.cs

namespace SpatialLab
{
    static class Program
    {
        private const double metersPerFoot = (2.54 * 12) / 100d;

        private static void DotSpatialTest1(int epsgCode, double x, double y, double radiusInMeters)
        {
            var sr1 = DotSpatial.Projections.ProjectionInfo.FromEpsgCode(epsgCode);
            DotSpatial.Data.IFeatureSet inputFeatures = new DotSpatial.Data.FeatureSet(DotSpatial.Topology.FeatureType.Point)
            {
                Projection = sr1,
            };

            string aeqd_proj_string = $"+proj=aeqd +lon_0={x} +lat_0={y} +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
            var sr2 = DotSpatial.Projections.ProjectionInfo.FromProj4String(aeqd_proj_string);

            DotSpatial.Topology.IBasicGeometry geom = new DotSpatial.Topology.Point(x, y);
            inputFeatures.AddFeature(geom);
            inputFeatures.Reproject(sr2);
            DotSpatial.Data.IFeatureSet outputFeatures = new DotSpatial.Data.FeatureSet(DotSpatial.Topology.FeatureType.Polygon)
            {
                Projection = inputFeatures.Projection
            };
            DotSpatial.Analysis.Buffer.AddBuffer(inputFeatures, radiusInMeters, outputFeatures);

            outputFeatures.Reproject(sr1);

            Log(outputFeatures.Extent.Width * metersPerFoot);
            Log(outputFeatures.Extent.Height * metersPerFoot);
        }

        static void OSGeoTest1(int epsgCode, double x, double y, double radiusInMeters)
        {
            // necessary to get gdal environment operational
            GdalConfiguration.ConfigureGdal();

            // define an input geometry
            var sr1 = new OSGeo.OSR.SpatialReference("");
            sr1.ImportFromEPSG(epsgCode);
            var g = OSGeo.OGR.Geometry.CreateFromWkt($"POINT ({x} {y})");
            g.AssignSpatialReference(sr1);

            // "define" a spatial reference
            var sr2 = new OSGeo.OSR.SpatialReference("");
            string aeqd_proj_string = $"+proj=aeqd +lon_0={g.GetX(0)} +lat_0={g.GetY(0)} +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
            sr2.ImportFromProj4(aeqd_proj_string);
            g.TransformTo(sr2);
            g = g.Buffer(radiusInMeters);
            g.TransformTo(sr1);

            OSGeo.OGR.Envelope env = new OSGeo.OGR.Envelope();
            g.GetEnvelope(env);
            double dx = env.MaxX - env.MinX;
            double dy = env.MaxY - env.MinY;
            Log(dx * metersPerFoot);
            Log(dy * metersPerFoot);
        }

        private static void Log(object message)
        {
            System.Console.WriteLine(message);
        }

        private static OSGeo.OGR.Geometry Buffer(this OSGeo.OGR.Geometry g, double meters)
        {
            var buffer = g.Buffer(meters, 8);
            buffer.AssignSpatialReference(g.GetSpatialReference());
            return buffer;
        }

        private static string AsWkt(this OSGeo.OGR.Geometry buffer)
        {
            string wkt;
            buffer.ExportToWkt(out wkt);
            return wkt;
        }

        static void Main(string[] args)
        {
            OSGeoTest1(x: 1571547, y: 1112049, epsgCode: 32033, radiusInMeters: 100);
            DotSpatialTest1(x: 1571547, y: 1112049, epsgCode: 32033, radiusInMeters: 100);
            System.Console.Read();
        }

    }
}