Seeking library for computing accurate geographic buffer around coordinate in EPSG:4326?

Question

In short, are there any C# or JavaScript API's that are capable of replacing https://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer/buffer?

I am trying to compute a geographically accurate buffer around a coordinate in EPSG:4326, but apart from using the ESRI Geometry Service I know of no API's to compute this buffer (maybe the Java version of GeoAPI can do it).

I asked a similar question but it was too specific to the .NET port of GeoAPI and has no answers.

Ideally a .NET library exists, but a JavaScript (Node or browser) solution would work, too.

I do not have the luxury of knowing which local coordinate system I should use, otherwise I could transform into that coordinate system and buffer using feet/meters depending on the unit-of-measure before transforming the entire buffer back to 4326.

I know of algorithms that will produce a point given a source location (lon/lat), initial bearing and distance (m) but I am trying to avoid implementing that myself.

I looked at turf-buffer but the code makes me think this is just an approximation optimized for the equator. When units are meters it simply divides the radius by 111120, note that 360 degrees * 111120 meters = 40003 km (the circumference of the earth).

As an alternative solution, perhaps there is a way of choosing a reasonable local coordinate system given a coordinate on 4326?

Answer 1

C#.Net Solution using GDAL/OGR bindings..

If you're willing to add the GDAL/OGR C# bindings to your .Net project (see second heading), then you can do what you want within a .Net environment. It's the same thing I did above in my prior answer using JavaScript, but 100% in C# here, and only one library— the GDAL C# bindings.

One (major) improvement of this implementation over my JavaScript example, is this version takes Well Known Text geometries, which can be point, line, polygon, multipart, etc., and buffers them without your having to additionally explode the individual coordinates then reassemble the geometries. So this is a huge improvement.

Like my JavaScript approach above, though, this route still applies the concept of using an Azimuthal Equidistant Projection centered on the geometry, then buffers on the localized grid, then transforms the geometry back to geographic coordinates. The projection is re-centered for every geometry submitted.

I drew from this GDAL C# example as a reference. Other GDAL C# examples are here.

using System;
using OSGeo.OSR;
using OSGeo.OGR;

namespace SharpGEO
{
    class GeoUTIL
    {
        public static string BufferGeom(string wkt)
        {
            string projected_geom_WKT = "";

            // If GDAL > 3.0..
            // Unexpected coordinate reversal in Well Known Text!
            // Explaination: https://github.com/OSGeo/gdal/issues/1546
            // Solution is.. 
            // SpatialReference.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);

            SpatialReference sr1 = new SpatialReference("");
            sr1.ImportFromEPSG(4326);
            sr1.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);

            Geometry g1 = Geometry.CreateFromWkt(wkt);
            g1.AssignSpatialReference(sr1);

            // Create localized Azimuthal Equidistant Projection centered on the point you want to buffer..
            // Example:
            // "+proj=aeqd +lon_0=-80.9957 +lat_0=33.97105 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=ft +no_defs";

            // This block uniquely adapts the Azimuthal Equidistant Projection for every geometry submitted.
            // We get the centroid in case a linestring or a polygon is submitted.
            // This keeps the projection centered even for non-point shapes.
            Geometry proj_center_geom = g1.Centroid();
            string aeqd_proj_string = "+proj=aeqd +lon_0=" + proj_center_geom.GetX(0).ToString() +
                                                " +lat_0=" + proj_center_geom.GetY(0).ToString() +
                                                " +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=ft +no_defs";

            // Transform to localized projection..
            SpatialReference sr2 = new SpatialReference("");
            sr2.ImportFromProj4(aeqd_proj_string);
            sr2.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);
            g1.TransformTo(sr2);

            // Buffer..
            // https://gdal.org/python/osgeo.ogr.Geometry-class.html#Buffer
            // 50 = distance in projection units (feet, in this case), and 30 = "quadsecs"
            // "QuadSecs" are the number of segments used to represent 90-degree arc of circle..
            Geometry g2 = g1.Buffer(50, 30);
            g2.AssignSpatialReference(sr2); // We have to tell OGR the projection of the new geom

            string g2WKT = "";
            g2.ExportToWkt(out g2WKT);
            Console.WriteLine(g2WKT);

            // Transform back to WGS84/LatLng
            g2.TransformTo(sr1);
            string g2GeogWKT = "";
            g2.ExportToWkt(out g2GeogWKT);
            //Console.WriteLine(g2GeogWKT);

            projected_geom_WKT = g2GeogWKT;
            return projected_geom_WKT;
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            // Point you want to buffer..
            string latlngWKT = "POINT(-80.9957 33.97105)";
            string proj_geom_WKT = GeoUTIL.BufferGeom(latlngWKT);

            Console.WriteLine("-------------------------------");
            Console.WriteLine(proj_geom_WKT);
        }
    }
}

Importing GDAL/OGR C# bindings into your C#.Net project..

I had to import the following GDAL/OGR assemblies (.dll files) into my C# project to reference the bindings:

gdal_csharp.dll ogr_csharp.dll osr_csharp.dll

I found them where I initially unpacked files for my GDAL installation in the following directory:

\bin\gdal\csharp

I also found them here (probably where GDAL copied them during install):

C:\Program Files (x86)\GDAL\csharp

One last thing, I ran into some build issues at first, and I had to go to Project > Properties > Build (tab), and set "Platform Target" to "x86".

A comment on that ..I probably intentionally installed 32-bit GDAL for macro compatibility with a lot of other things. But for this purpose, it might make sense to download a specific GDAL, like a 64-bit build, for importing into your project.

Answer 2

Use GeographicLib (NETGeographicLib, or JavaScript, or other implementation).

You would then need to create a method that would:

• Use an ellipsoid of revolution definition, usually WGS84 for EPSG:4326
• Take a central point lat1, lat2, and buffer distance s12 in metres
• Compute a series of direct geodesic calculations by varying azi1 from 0 to 360 degrees, i.e. "draw a circle"
• From each direct geodesic calculation, take lat2 and lon2 for the exterior ring of the buffered geometry.

The round-off errors in the direct geodesic calculation is less than 15 nanometers.

Answer 3

You might consider 1: defining an Azimuthal Equidistant projection centered on the point you want to buffer, 2: reprojecting your coordinate from LatLong to the custom azimuthal projection, then 3: buffering the point using whatever library you prefer against the Azimuthal Equidistant coordinates, which preserves distance so should give you a clean buffer, then 4: iterate back over the new coordinates and convert them to the desired projection

In a browser environment, I created an example of this approach posted on JSFiddle, using the JavaScript Topology Suite and Proj4JS libraries. You can see it here:

http://jsfiddle.net/aLvp0gqm/3/

If you wanted to stay in .Net, you could take the approach discussed above (the essence of which is usage of a centered Azimuthal Equidistant projection) and adapt the guts of this solution, which uses the .Net port of the Proj.4 project to convert between coordinate systems.

Ideally you could do all this with one library but I don't have an adaptable solution ready to go from my bag of tricks. One option I'd consider, though, would be looking for a solution using the OGR .Net bindings. Honestly that would probably be a cleaner solution. Maybe I can return to this later and experiment with that angle.