Fun Stuff: Computing circular buffers in geographic coordinates
Finding all the objects within a certain distance from a point is surely a common GIS problem. The problem is normally solved using OGC “dwithin” filters or by computing a buffer and then finding all the intersecting objects.
Very often both of the approaches fail miserably in case the coordinate system is a geographic one, as common libraries, such as JTS and GEOS, are not able to handle the non planar nature of it. As far as “dwithin” is concerned rencent Oracle and PostGIS versions can manage the problem properly, but what to do if they cannot be used?
We had to solve this problem when computing data distribution statistics over raster data cells that are within a certain distance from a given point, and making for an accurate calculation regardless of how long the distance was.
To do that we created a new GeoServer WPS process, “gs:PointBuffers”, that can create a set of buffers given a point, a target SRS and a set of distances in meters.
In case the SRS denotes a geographic spatial reference system the GeoTools GeodeticCalculator is used to sample the set of points that are at the given distance, looping over a closed sets of azimuths to cover the entire shape.
Interested in seeing the results? I certainly was.
Let’s start with a set of small buffers at a medium latitude: 10, 30, 50 and 100 km buffers around a point located in northern Italy. Here is there result:
As you can see, drawing the result in plain WGS84 (plate carré for the conoisseurs) we get elliptical shapes. This should not come as a surprise if you consider that at 45° one degree of latitude spans 111km, whilst a degree of longitude spans only 78km (see the “Degree length” table at Wikipedia).
What if we pump up the distance significantly? Let’s try with 100, 500, 1000, 2000 and 3000km instead. Here is the result:
See the funny shape we get? This is the effect of the size of one degree of longitude shrinking as we move towards north.
It is also a good indicator of how deformed the now common WGS84 maps, often published on the web, are.
If you want to see the same data in a common projection, let’s have a look at the same map in EPSG:3857 (aka the Google projection):
Somewhat better, even if the Mercator tendency to inflate areas at high latitudes is well evident.
Well, this is it. The gs:PointBuffers is soon going to land in GeoServer for your testing pleasure.
We’d very much like to tackle the same problem against lines and polygons as well. Interested? Let us know!
The GeoSolutions team