How do I create a circle with latitude, longitude and radius with GeoTools?

Question

Right now I have:

Polygon circle = geometryBuilder.circle(
myLong,
myLat, 
radiusInMeters, 10);

And it creates (with lat=28.456306, long=-16.292034 and radius=500) a nonsense polygon with huge latitudes and longitudes, such as:

POLYGON ((483.678055 28.482505000000003, 388.1865521874737 -265.4101211462366, 138.1865521874737 -447.04575314757676, -170.8304421874737 -447.0457531475768, -420.8304421874737 -265.41012114623663, -516.321945 28.482504999999943, -420.83044218747375 322.3751311462365, -170.8304421874738 504.01076314757677, 138.18655218747358 504.0107631475768, 388.18655218747364 322.3751311462367, 483.678055 28.482505000000003))

I expected to have ten pairs of coordinates with lat's and long's nearby the center point I supplied.

Any help would be more than helpful. Thanks in advance!

EDIT

In addition to @iant 's answer, I had to create a Point as a Feature

//build the type
SimpleFeatureType TYPE = null;
try {
    TYPE = DataUtilities.createType("", "Location", "locations:Point:srid=4326," + "id:Integer" // a
            // number
            // attribute
            );
} catch (Exception e1) {
    // TODO Auto-generated catch block
    e1.printStackTrace();
}

SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
com.vividsolutions.jts.geom.Point point = geometryFactory.createPoint(
        new Coordinate(
                currentDevicePosition.getLongitude(), 
                currentDevicePosition.getLatitude()
                )
        );
featureBuilder.add(point);
SimpleFeature feature = featureBuilder.buildFeature( "fid.1" ); // build the 1st feature

as explained in iant's Gist here: https://gitlab.com/snippets/17558 and here: http://docs.geotools.org/, oh, and also I was missing a dependency as explained here SchemaException in java

Samalla

Answer 1

There are two solutions to this:

1. Convert the radius in meters to degrees and treat the problem as a planar problem

2. Convert the lat/lon point to meters, calculate the circle in a locally planar projection and reproject back to lat/lon.

For 1 you could do something like which will be fine for small radii near the equator:

GeodeticCalculator calc = new  GeodeticCalculator(DefaultGeographicCRS.WGS84);
  calc.setStartingGeographicPoint(point.getX(), point.getY());
  calc.setDirection(0.0, 10000);
  Point2D p2 = calc.getDestinationGeographicPoint();
  calc.setDirection(90.0, 10000);
  Point2D p3 = calc.getDestinationGeographicPoint();

  double dy = p2.getY() - point.getY();
  double dx = p3.getX() - point.getX();
  double distance = (dy + dx) / 2.0;
  Polygon p1 = (Polygon) point.buffer(distance);

I'll show some code for the second as it is more general (i.e. it works better and for a better range of radii).

First you need to find a local projection, GeoTools provides a "pseudo" projection AUTO42001,x,y which is a UTM projection centred at X,Y:

public SimpleFeature bufferFeature(SimpleFeature feature, Measure<Double, Length> distance) {
    // extract the geometry
    GeometryAttribute gProp = feature.getDefaultGeometryProperty();
    CoordinateReferenceSystem origCRS = gProp.getDescriptor().getCoordinateReferenceSystem();

    Geometry geom = (Geometry) feature.getDefaultGeometry();
    Geometry pGeom = geom;
    MathTransform toTransform, fromTransform = null;
    // reproject the geometry to a local projection
    if (!(origCRS instanceof ProjectedCRS)) {

      double x = geom.getCoordinate().x;
      double y = geom.getCoordinate().y;

      String code = "AUTO:42001," + x + "," + y;
      // System.out.println(code);
      CoordinateReferenceSystem auto;
      try {
        auto = CRS.decode(code);
        toTransform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, auto);
        fromTransform = CRS.findMathTransform(auto, DefaultGeographicCRS.WGS84);
        pGeom = JTS.transform(geom, toTransform);

      } catch (MismatchedDimensionException | TransformException | FactoryException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }

    }

So now pGeom is our point in metres. Buffering it is now easy

  Geometry out = bufferFeature(pGeom, distance.doubleValue(SI.METER)); 

then we project back to WGS84 (lat/lon) using the reverse transform we looked up earlier:

  retGeom = JTS.transform(out, fromTransform); 

There is then a little messing around to change the feature type to reflect the fact we are returning a polygon instead of a Point. The full code is in this gist.

When I run it I get the following output:

POINT (10.840378413128576 3.4152050343701745)
POLYGON ((10.84937634426605 3.4151876838951822, 10.849200076653755 3.413423962919184, 10.84868480171117 3.4117286878605766, 10.847850322146979 3.4101670058279794, 10.846728706726902 3.4087989300555464, 10.845363057862208 3.407677033830687, 10.843805855306746 3.406844430298736, 10.84211693959797 3.406333115754347, 10.840361212705258 3.4061627400701946, 10.838606144204721 3.4063398515107184, 10.836919178768184 3.4068576449605277, 10.835365144548726 3.4076962232621035, 10.834003762019957 3.408823361646906, 10.832887348980522 3.410195745914279, 10.832058809914859 3.411760636805914, 10.831549986992338 3.4134578966399034, 10.831380436105858 3.4152223003379722, 10.831556675029052 3.416986042039048, 10.832071932633442 3.4186813409639054, 10.832906408849936 3.4202430463705085, 10.834028035422469 3.4216111414662183, 10.835393708241908 3.422733050021835, 10.836950943907517 3.4235656570147763, 10.838639896841123 3.424076965623486, 10.840395659406198 3.4242473268789406, 10.842150756595839 3.4240701947133396, 10.843837739370569 3.4235523773972796, 10.845391776937724 3.4227137757216988, 10.846753148314034 3.4215866180136185, 10.847869537398722 3.4202142214154887, 10.848698043354238 3.4186493270628633, 10.849206829051935 3.4169520731645546, 10.84937634426605 3.4151876838951822))

Fung

Having copied the whole function and the auxiliar buffer from your gist, how do I call it having a polygon (or a GeoJSON definition geojson.org/geojson-spec.html#polygon of it) instead of a feature? – Samalla

I think but want to confirm that by changing to projection and then back again, will account for the differences in meters / latitude at different longitudes? – Battles

Also, #1 would be more performant? What if you had to do this 1M times? – Battles

@Samalla you could just skip the part before geom is extracted from the feature – Fung

@Battles - yes this handles the projection issues for you. #1 would be slightly faster (you still need to adjust the radius for each point) but less accurate for any realistically sized buffer – Fung

thanks @iant, anyway, could you ellaborate more on #1? I'm fairly new to Geotools, but anyway, how I change "a radius in meters to degrees"? AFAIK, a should be a pair of points, or a magnitude, and I don't know how to convert it to a "double radius" – Samalla

what are the packages used? i can't find the point.buffer() method??? also , how did you get the 'distance' in #1 – Tattoo

@JerylCook - it should be bufferFeature - I've fixed it. – Fung

the answer is a bit confusing..can you clean it up with just the solution at the top? ur answer 'Polygon p1 = (Polygon) point.buffer(distance);' , this .buffer is not implemented..and you said use bufferFeature doesn't take 'distance' as a parameter. – Tattoo

It has been a while since the answer, so I would have a small question. Is there a better way than approach nr. 2 ? maybe there is a library (that I am not aware of) which already provides us with the whole mechanism ? – Swollen

spatial4j (github.com/locationtech/spatial4j) has nice support for circle and you convert it to jts geometry – Perambulator

Answer 2

    double latitude = 40.689234d;
    double longitude = -74.044598d;
    double diameterInMeters = 2000d; //2km

    GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
    shapeFactory.setNumPoints(64); // adjustable
    shapeFactory.setCentre(new Coordinate(latitude, longitude));
    // Length in meters of 1° of latitude = always 111.32 km
    shapeFactory.setWidth(diameterInMeters/111320d);
    // Length in meters of 1° of longitude = 40075 km * cos( latitude ) / 360
    shapeFactory.setHeight(diameterInMeters / (40075000 * Math.cos(Math.toRadians(latitude)) / 360));

    Polygon circle = shapeFactory.createEllipse();

Bergman

Is longitude and latitude mixed up in the centre coordinate? – Nightrider

Yeah I think it is – Nine

If the shape was a rotated ellipse say to 70 degrees clockwise, will this also be correct in that case? – Telefilm

Answer 3

GeometricShapeFactory shapeFactory = new GeometricShapeFactory(); shapeFactory.setNumPoints(64); adjustableshapeFactory.setCentre(new Coordinate(latitude, longitude)); // Length in meters of 1°

Haskins

Please explain a bit your answer. – Respire