geodetic buffers with geotools

sometimes you need a sphere or even better an ellipsoid…

the above images shows several point buffers all with a radius of 500 km… web mercator at it’s best…

Following two quick hacks for a Maven/geotools environment. For more details on setting up geotools with maven look at this tutorial…

simple java application producing GeoJSON polygon using geotools JSON encoder

this example produces a point buffer, but you can use the built in line-buffer as well. A sample linestring (line) is included… Note that the linebuffer- implementation is just a draft… far from correct, but it delivers robust results for straight lines.

import java.awt.geom.Point2D;
import java.io.IOException;
import java.io.StringWriter;
import java.util.ArrayList;
import java.util.Collections;

import javax.xml.transform.TransformerException;

import org.geotools.feature.DefaultFeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geojson.feature.FeatureJSON;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.GeodeticCalculator;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.xml.sax.SAXException;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.LinearRing;
import com.vividsolutions.jts.geom.Polygon;

public class TestGeodetic {

 public static Geometry getGeodeticPointBuf(Double lon, Double lat, int rad)
 {
 //calculate on WGS84 ellipsoid
 GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
 calc.setStartingGeographicPoint(lon, lat);
 GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null );
 ArrayList<Coordinate> coords = new ArrayList<Coordinate>();

 //create 50 points with polar coordinates (angle+distance) and add to ArrayList - TODO increase point amount when radius increases...
 for (int i = 0; i < 50; i++)
 {
 int angle = (i * 360 / 50);
 if (angle>=180){angle = angle-360;}
 calc.setDirection(angle, rad);
 Point2D dest = calc.getDestinationGeographicPoint();
 coords.add(new Coordinate(dest.getX(),dest.getY(),0));

 }
 coords.add(coords.get(0)); //close ring

 Coordinate[] coords1 = new Coordinate[coords.size()];
 coords1 = coords.toArray(coords1);

 LinearRing ring = geometryFactory.createLinearRing( coords1 );
 LinearRing holes[] = null; // use LinearRing[] to represent holes
 Polygon polygon = geometryFactory.createPolygon(ring, holes );

 return polygon;
 }

 public static Geometry getGeodeticLineBuf(Geometry inline, int dist)
 {
 GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null );
 GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
 Coordinate[] subsat_points = inline.getCoordinates();

 ArrayList<Coordinate> hull_right = new ArrayList<Coordinate>();
 ArrayList<Coordinate> hull_left = new ArrayList<Coordinate>();
 for (int i = 0; i < subsat_points.length-1; i++)
 {

 calc.setStartingGeographicPoint(subsat_points[i].x,subsat_points[i].y);
 calc.setDestinationGeographicPoint(subsat_points[i+1].x,subsat_points[i+1].y);
 double angle = calc.getAzimuth();

 if (angle<=-90 ){angle = angle+180;}
 if (angle>=90 ){angle = angle-180;}
 calc.setStartingGeographicPoint(subsat_points[i].x,subsat_points[i].y);
 calc.setDirection(angle+90,dist);
 Point2D dest1 = calc.getDestinationGeographicPoint();
 hull_right.add(new Coordinate(dest1.getX(),dest1.getY(),0));
 calc.setDirection(angle-90, dist);
 Point2D dest2 = calc.getDestinationGeographicPoint();
 hull_left.add(new Coordinate(dest2.getX(),dest2.getY(),0));

 }
 Collections.reverse(hull_left);

 ArrayList<Coordinate> hull = new ArrayList<Coordinate>();
 hull.addAll(hull_right);
 hull.addAll(hull_left);
 hull.add(hull_right.get(0)); //close ring
 Coordinate[] hull1 = new Coordinate[hull.size()];
 hull1 = hull.toArray(hull1);

 LinearRing ring = geometryFactory.createLinearRing( hull1 );

 LinearRing holes[] = null; // use LinearRing[] to represent holes
 Polygon polygon = geometryFactory.createPolygon(ring, holes );

 return polygon;

 }

public static void main(String[] args) throws IOException, SAXException, TransformerException
{

 //create test-line... for linebuffers
 GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory( null );

Coordinate[] coords =
 new Coordinate[] {new Coordinate(-71, 52.6), new Coordinate(-56, 66), new Coordinate(-49, 73),new Coordinate(-40, 75), new Coordinate(-35, 85) };
 @SuppressWarnings("unused")
 LineString line = geometryFactory.createLineString(coords);

 SimpleFeatureTypeBuilder b = new SimpleFeatureTypeBuilder();

//set the name
 b.setName( "geodetic-buffer" );

//add some properties - just a template
 //b.add( "a", String.class );
 //b.add( "b", String.class );
 //b.add( "c", String.class );
 //b.add( "d", String.class );

 //add a geometry property
 b.setCRS( DefaultGeographicCRS.WGS84 ); // set crs first
 b.add( "geom", Polygon.class ); // then add geometry

//build the type
 final SimpleFeatureType FLAG = b.buildFeatureType();
 //
 SimpleFeatureBuilder builder = new SimpleFeatureBuilder(FLAG);

DefaultFeatureCollection poly_out = new DefaultFeatureCollection("internal",FLAG);
 //Line buffer
 //builder.add(getGeodeticLineBuf(line, 200000));

 //point buffer with 500000 meters
 builder.add(getGeodeticPointBuf(-50.3333, 66.5, 500000));

 SimpleFeature feature_out = builder.buildFeature(null);
 poly_out.add(feature_out);

 FeatureJSON fjson = new FeatureJSON();
 fjson.isEncodeFeatureCollectionCRS();
 StringWriter writer = new StringWriter();
 fjson.writeFeatureCollection(poly_out, writer);

 System.out.println(writer.toString().replace("\"properties\":", "\"properties\":{")); //tweak to produce valid geoJSON

}

}

geodetic point buffer as servlet without geotools JSON encoder (faster)

this servlet produces some GeoJSON for usage in a webapp e.g. ol3… Your request could look like http://yourdomain.com/GeodeticCircle?&lon=50.05257&lat=%20-27.91024&dist=500

response:

{
 "type" : "Polygon",
 "coordinates" : [[[50.05257, -23.39691514099216], [50.81729625396115, -23.450400161188764], [51.565142999703205, -23.60966376988598], [52.27950095965239, -23.87115272264491], [52.94430066379978, -24.22901487736714], [53.54428141982014, -24.675202669272224], [54.065260902584455, -25.19961551117886], [54.49440783160426, -25.790279456676938], [54.82052095616245, -26.433562922038107], [55.03431721008193, -27.114428081084437], [55.12872984685988, -27.816718485923438], [55.099213192896336, -28.523484166216896], [54.94404424451737, -29.217345478220302], [54.664603059566836, -29.880895852665123], [54.265604801971044, -30.497140951777368], [53.75524819424773, -31.04996743309319], [53.1452404735973, -31.524628742137907], [52.45066039974901, -31.908228813772958], [51.68963051111196, -32.19017844843427], [50.882788164856144, -32.362594977801486], [50.05257, -32.42061515391178], [49.222351835143854, -32.362594977801486], [48.415509488888034, -32.19017844843427], [47.654479600250994, -31.908228813772958], [46.95989952640271, -31.524628742137907], [46.34989180575228, -31.04996743309319], [45.839535198028955, -30.497140951777368], [45.44053694043317, -29.880895852665123], [45.161095755482634, -29.217345478220302], [45.005926807103656, -28.523484166216896], [44.97641015314012, -27.816718485923438], [45.070822789918076, -27.114428081084437], [45.28461904383756, -26.433562922038107], [45.61073216839574, -25.790279456676938], [46.03987909741554, -25.19961551117886], [46.56085858017986, -24.675202669272224], [47.16083933620022, -24.22901487736714], [47.82563904034762, -23.87115272264491], [48.53999700029679, -23.60966376988598], [49.28784374603885, -23.450400161188764], [50.05257, -23.39691514099216]]],
 "crs" : {
 "type" : "name",
 "properties" : {
 "name" : "urn:ogc:def:crs:OGC:1.3:CRS84"
 }
 }
}

import java.awt.geom.Point2D;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.ArrayList;

import javax.servlet.ServletException;
import javax.servlet.http.HttpServlet;
import javax.servlet.http.HttpServletRequest;
import javax.servlet.http.HttpServletResponse;
import org.geotools.referencing.GeodeticCalculator;
import org.geotools.referencing.crs.DefaultGeographicCRS;
public class GeodeticCircle extends HttpServlet {
 private static final long serialVersionUID = 1L;
 public GeodeticCircle() {super();}

 

protected void doGet(HttpServletRequest request, HttpServletResponse response) throws ServletException, IOException {
 response.setContentType("text/plain");
 response.setCharacterEncoding("UTF-8");
 PrintWriter printWriter = response.getWriter();

 Double lon = Double.parseDouble(request.getParameter("lon"));
 Double lat = Double.parseDouble(request.getParameter("lat"));
 int distance = Integer.parseInt(request.getParameter("dist"))*1000;

 GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
 calc.setStartingGeographicPoint(lon, lat);

 ArrayList<String> coords = new ArrayList<String>();

 for (int i = 0; i < 40; i++)
 {
 int angle = (i * 360 / 40);
 if (angle>=180){angle = angle-360;}
 calc.setDirection(angle /* azimuth */, distance/* distance */);
 Point2D dest = calc.getDestinationGeographicPoint();
 coords.add("["+Double.toString(dest.getX())+","+Double.toString(dest.getY())+"]");
 }
 coords.add(coords.get(0)); //close ring

 printWriter.println("{\"type\":\"Polygon\",\"coordinates\":["+coords+"],\"crs\":{\"type\":\"name\",\"properties\":{\"name\":\"urn:ogc:def:crs:OGC:1.3:CRS84\"}}}");
 }

protected void doPost(HttpServletRequest request, HttpServletResponse response) throws ServletException, IOException {
 doGet(request, response);
 }

}