如何使用GeoTools创建带有纬度、经度和半径的圆形？

Question

现在我有：

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

它使用纬度为28.456306、经度为-16.292034和半径为500创建了一个无意义的多边形，包含巨大的纬度和经度，例如：

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

我希望我提供的中心点附近有十对纬度和经度的坐标。

任何帮助都将非常有用。提前感谢！

编辑

除了 @iant 的答案外，我还需要创建一个作为要素的点。

//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

正如iant在这里的Gist中所解释的那样：https://gitlab.com/snippets/17558和http://docs.geotools.org/，哦，还有我错过了一个依赖项，就像这里所解释的那样：Java中的SchemaException

- jonayreyes

Answer 1

有两种解决方案：

1. 将半径从米转换为度数，将问题视为平面问题。

2. 将纬度/经度点转换为米，在本地平面投影中计算圆圈，然后重新投影回纬度/经度。

对于第一种方法，你可以做类似这样的事情，适用于赤道附近小半径的情况:

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

因为第二种方法更通用（即它的适用范围更广，效果也更好），所以我将展示一些代码。

首先，您需要找到一个本地投影，GeoTools提供了一个“伪”投影AUTO42001,x,y，它是以X,Y为中心的UTM投影：

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();
      }

    }

所以现在pGeom是我们的点，单位为米。对它进行缓冲现在很容易。

Translated content:

因此，现在pGeom是以米为单位的点。 对其进行缓冲现在很容易。

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

然后，使用之前查找的反向变换将其投影回WGS84（纬度/经度）：

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

  retGeom = JTS.transform(out, fromTransform); 

然后需要进行一些调整，以更改要返回的特征类型以反映我们返回的是多边形而不是点。完整代码在此代码片段中。

运行代码后，将得到以下输出：

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

- Ian Turton

---

@jonayreyes 你可以直接跳过从要素中提取 geom 之前的部分。 - Ian Turton

谢谢@iant，不过，你能详细解释一下第一点吗？我对Geotools还比较新，但是，我该如何将“以米为单位的半径转换为度数”？据我所知，a应该是一对点或一个大小，我不知道如何将其转换为“双倍半径”。 - jonayreyes

用了哪些包？我找不到 point.buffer() 方法？？？另外，你是怎么得到 #1 中的“distance”的？ - Jeryl Cook

答案有点混乱，你能否只给出顶部的解决方案？你的答案 'Polygon p1 = (Polygon) point.buffer(distance);'，这个 .buffer 方法没有实现，而且你说使用 bufferFeature 但它不接受 'distance' 参数。 - Jeryl Cook

spatial4j（https://github.com/locationtech/spatial4j）对圆形有很好的支持，您可以将其转换为JTS几何图形。 - Piotr

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

- Jared

---

中心坐标中经度和纬度是否混淆了？ - ThomasRS

是的，我认为是的。 - Sparm

Answer 3

GeometricShapeFactory shapeFactory = new GeometricShapeFactory(); shapeFactory.setNumPoints(64); shapeFactory.setCentre(new Coordinate(latitude, longitude)); // 1°的长度（米）

- Thirty Axis

---

请稍微解释一下您的答案。 - Kaveh