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

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?

## 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 sr1 = new SpatialReference("");
sr1.ImportFromEPSG(4326);

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

• I'm sorry to be so dense...so you are creating a reference system centered near the feature of interest that way you know 1 degree north is about the same as 1 degree west, so it is safe to assume there is a simple conversion from meters to degrees (360 degrees/40003 meters). Once you have your circle built you transform it back to the desired projection and you get an geometry such that every point on the border is equidistant to the centroid. Is that a decent summary?
– ca0v
Jul 7, 2021 at 13:54
• (1 of 2): One thing that doesn't feel quite right about your summary is this bit: ..1 degree north is about the same as 1 degree west, so it is safe to assume there is a simple conversion from meters to degrees. Once projected, you're using cartesian units. (Aside: you'll want to be mindful if you want your units to be meters or feet. I used feet, above, as indicated by `+units=ft` in the projection definition. If you want meters, use `+units=m` ..that's very important.) Like all projections, some traits are preserved, and other traits are distorted. But.. Jul 7, 2021 at 14:14
• (2 of 2): In the case of this projection, distance is preserved measuring away from the center of the proj. So that's perfect for buffering away from the center, especially for buffering a point. I expect it's also good for buffering a shape, providing it's centered in the same way. If you look at the Tissot indicatrix for this, areal distortion is extraordinary near the perimeter. But close to the center it's very good. Jul 7, 2021 at 14:24
• Ha! Trivial diversion, but I just noticed you're in Greenville SC. I'm in Columbia SC. Small world :))) Jul 7, 2021 at 14:27
• It was all there in front if me, I just needed to hear it twice. I now have a working GDAL sample. If you make it to Greenville I definitely owe you a beer/coffee/mountain bike guided tour.
– ca0v
Jul 7, 2021 at 14:31

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.

• I noticed NETGeographicLib has no methods that begin with "b", including "buffer" but your "draw a circle" solution would work. I will give this a try.
– ca0v
Jul 1, 2021 at 22:35

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.

• How would I know which local projection you use? How did you determine "ELROBIS:101" was suitable for coordinate you selected? Will that coordinate system work for any coordinate in the world? I do not know how to select a suitable projection given a coordinate.
– ca0v
Jul 1, 2021 at 22:30
• Ha! Ok so the "ELROBIS:101" name was just a joke. It's actually a custom projection, but it's an easy one to implement. If you check the comment right above that in the fiddle example, it should mention that you can name it whatever you want (abracadabra, ca0v:1, etc), and you can easily re-center it by just changing the latitude longitude values to the values of the actual point you're buffering. Then, all the distance calculations away from that point will always be true. Hopefully that makes sense. Jul 1, 2021 at 22:54
• To explain about more, if you were to buffer 1000 different points using this approach, you'd ultimately redefine the projection 1000 different times, each time centered on the point that you were actually buffering. That might sound like it would be slow, but I don't think it would be. It's essentially just changing the projection parameters. Jul 1, 2021 at 22:56
• One more thought, if you wanted to experiment with the OGR .Net bindings I mentioned in my last paragraph, I found an old answer where I provided a similar example using C#/.Net, it might be a little more involved, but you can look it over here: gis.stackexchange.com/a/61574 Jul 1, 2021 at 22:59
• I am concerned that your leaflet example produces what looks to be a circle against a 3857 projection. If the distance measurements were accurate wouldn't I expect to see an oval? In fact I do see an oval when I use the 0,0 projection. It seems the proj4.buffer is doing exactly what I want without using a custom projection.
– ca0v
Jul 6, 2021 at 14:16