基于Java和PostGIS的AOI面数据球面面积计算实践

前言

一、计算方法简介

二、球面面积计算

1、AOI数据转Polygon

2、Geotools面积计算

3、GeographicLib面积计算

4、PostGIS面积计算

三、结果分析

1、不同算法结果对比

2、与互联网AOI对比

3、与天地图测面对比

四、总结

前言

在现代地理信息系统（GIS）中，AOI（Area of Interest）面数据的处理和分析是一个核心任务。AOI面数据通常用于表示地理区域，如城市、国家、自然保护区等。准确计算这些区域的面积对于城市规划、环境监测、资源管理等领域至关重要。然而，由于地球的曲率，传统的平面面积计算方法在处理大范围地理数据时会产生显著误差。在实际应用中，AOI面数据的球面面积计算可以应用于多种场景。例如，在城市规划中，准确计算城市区域的面积有助于评估城市扩展的潜力；在环境监测中，计算自然保护区的面积有助于评估生态系统的健康状况；在资源管理中，计算农田的面积有助于评估农业生产的潜力。因此，基于Java和PostGIS的AOI面数据球面面积计算具有广泛的应用前景。因此，球面面积计算成为解决这一问题的关键。

在Java中实现球面面积计算时，可以使用JTS Topology Suite（JTS）库来处理空间数据。JTS提供了丰富的几何操作函数，支持点、线、面等空间数据类型的处理。通过JTS，可以方便地读取AOI面数据，并将其转换为球面坐标系。然后，可以应用球面三角剖分方法，计算每个球面三角形的面积，并累加得到总面积。PostGIS的ST_Area函数可以直接计算球面面积，但为了验证自定义算法的准确性，可以将Java计算的结果与PostGIS的结果进行对比。通过对比，可以发现两者之间的差异，并进一步优化自定义算法。

本文详细讲解如何基于Java语言，分别使用Geotools、GeographicLib和PostGIS空间函数来分别计算对应的AOI数据的球面面积，为了验证计算结果。我们将采用天地图标绘、互联网查询公开数据的方式来验证空间数据的范围面积。天地图标绘由于采样的点趋向于简单，因此可能采集的精度可能有所损失。而使用互联网的数据则由于可能采集的面范围不一致，存在一定的误差。本例将说明三种不同的方式与实际结果的误差值。为大家对面状的球面面积计算有一个参考。

一、计算方法简介

在实际应用中，AOI面数据通常以矢量数据的形式存储，如Shapefile、GeoJSON等。这些数据包含了地理区域的边界信息，通常由一系列经纬度坐标点组成。为了计算球面面积，需要将这些坐标点转换为球面坐标系，并应用球面几何公式进行计算。PostGIS提供了ST_Area函数，可以直接计算球面面积，但为了更深入地理解计算过程，可以通过Java实现自定义的球面面积计算算法。球面面积计算的核心在于将地理坐标转换为球面坐标系，并应用球面三角剖分方法。球面三角剖分将AOI面数据划分为多个球面三角形，然后计算每个三角形的面积并累加得到总面积。球面三角形的面积计算可以使用L'Huilier公式或球面余弦定理。这些公式考虑了地球的曲率，能够提供更准确的计算结果。

为了确保球面面积计算的准确性，可以采用多种方法进行验证。例如，可以使用已知面积的地理区域进行测试，比较计算结果与已知面积的差异。此外，还可以使用不同的算法进行对比，如平面面积计算与球面面积计算的对比，以评估球面面积计算的必要性。需要说明的是，在面向GIS的应用当中，对于空间坐标系可以分为地理坐标系和投影坐标系。地理坐标系一般用经纬度标识，它的单位是度，因此直接进行面积计算得出的单位是平方度，是现实中不太好理解的。而投影坐标系是米制单位，通常符合实际的计算结果，如果数据已经是投影坐标，可以直接进行计算就可以得到平方米的面积。

二、球面面积计算

本节将重点介绍在Java和PostGIS空间数据库中来分别进行球面面积计算的方法。计算之前，我们准备了6组对比的数据，这6组数据来自于高德地图中获取的AOI数据。后续我们将会对这六组数据来进行球面面积求解，在后续的对比中也会对这几组数据进行对比。

1、AOI数据转Polygon

在高德地图中，AOI数据是一个收尾坐标连接起来的坐标字符串。因此我们需要对这串坐标字符进行解析，转换成我们的空间面数据，后续针对这类数据来进行面积求解。因此首先我们要将字符串的数据转换成Polygon，这里我们以Geotools为例，讲解如何把字符串数组转换成Polygon，核心代码如下：

/**
* - 将AOI字符串转换成Polygon对象
* @return
*/
public static Polygon convertAoi2Polygon(String aoistr) {
	String [] AOI_Str_Array = aoistr.split(";");
	// 获取GeometryFactory实例
    GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
    Coordinate[] coords = {};
    if( AOI_Str_Array.length > 0) {
        coords = new Coordinate[AOI_Str_Array.length];
    }
    //处理坐标
    for (int i = 0; i < AOI_Str_Array.length; i++) {
        String loc = AOI_Str_Array[i];
        String [] latlon = loc.split(",");
		double lng = Double.parseDouble(latlon[0]);
		double lat = Double.parseDouble(latlon[1]);
		//将高德坐标转换成WGS84坐标
		double [] gcj284 = CoordinateTransformUtil.gcj02towgs84(lng, lat);
		coords[i] = new Coordinate(gcj284[0], gcj284[1]);
	}
    LinearRing shell = geometryFactory.createLinearRing(coords);
    Polygon polygon = geometryFactory.createPolygon(shell, null);
    polygon.setSRID(4326);//设置4326的坐标
    return polygon;
}

2、Geotools面积计算

这里介绍使用Geotools来实现简化的球面面积近似计算，当然这里的地球半径等使用的一个固定常数，6371008.8。首先需要在Maven中引入以下依赖：

<dependency>
	<groupId>org.geotools</groupId>
	<artifactId>gt-referencing</artifactId>
	<version>28.2</version>
</dependency>

然后进行球面面积的近似计算，核心方法如下：

// 简化的球面面积近似计算（实际项目需优化）
private static double approximateSphericalArea(Coordinate p1, Coordinate p2) {
  // 此处为示例逻辑，真实计算需使用积分或地理公式
  // 例如：使用 Haversine 公式计算梯形面积
  double R = 6371008.8; // 地球平均半径（米）
  double dLon = Math.toRadians(p2.x - p1.x);
  double lat1 = Math.toRadians(p1.y);
  double lat2 = Math.toRadians(p2.y);
  return 0.5 * R * R * dLon * (Math.sin(lat2) + Math.sin(lat1));
}

最后定义一个使用 GeodeticCalculator 逐边计算多边形面积的方法，在这个方法的实现中，我们使用累加的方法来进行计算，核心方法如下：

/**
 * -使用 GeodeticCalculator 逐边计算多边形面积
 */
public static double calculateGeodeticArea(Geometry geometry, CoordinateReferenceSystem crs) throws Exception {
   if (!(geometry instanceof Polygon)) {
      throw new IllegalArgumentException("仅支持 Polygon 类型");
   }
   Polygon polygon = (Polygon) geometry;
   Coordinate[] coords = polygon.getExteriorRing().getCoordinates();

   GeodeticCalculator calculator = new GeodeticCalculator(crs);
   double totalArea = 0.0;

   // 使用球面三角形累加算法（类似 JTS 原有逻辑）
   for (int i = 0; i < coords.length - 1; i++) {
       Coordinate p1 = coords[i];
       Coordinate p2 = coords[i + 1];
       calculator.setStartingGeographicPoint(p1.x, p1.y);
       calculator.setDestinationGeographicPoint(p2.x, p2.y);
       // 累加每个边的贡献面积（需自行实现积分逻辑，此处简化）
       // 实际应参考 JTS 原有算法或使用第三方库（如 GeographicLib）
       totalArea += approximateSphericalArea(p1, p2);
    }
   return Math.abs(totalArea); // 返回绝对值
}

3、GeographicLib面积计算

如果觉得使用Geotools来计算面积比较繁琐，接下来我们使用一个简化的库，利用这个库可以更便捷的对一个面数据的面积进行求解。这个库就是：GeographicLib-Java。首先依然需要的Pom.xml中引用以下依赖：

<dependency>
	<groupId>net.sf.geographiclib</groupId>
	<artifactId>GeographicLib-Java</artifactId>
	<version>1.52</version>
</dependency>

计算球面面积的核心方法更加简单，如下所示：

public static double calculateExactGeodeticArea(Polygon polygon) {
	Geodesic geodesic = Geodesic.WGS84;
	PolygonArea areaCalc = new PolygonArea(geodesic, false);
	for (Coordinate coord : polygon.getExteriorRing().getCoordinates()) {
	     areaCalc.AddPoint(coord.y, coord.x); // 注意经纬度顺序
	}
	PolygonResult result = areaCalc.Compute();
	return Math.abs(result.area); // 单位：平方米
}

4、PostGIS面积计算

在空间数据库中进行球面面积的计算比较简单，为了求解平方米的单位，我们可以将数据转换成geograghy的类型来进行球面的计算，在使用空间面积函数st_area即可实现面积的求解，这里以某小学的空间面数据为例，讲解如何在PostGIS数据库中实现球面面积求解。具体可以体现为以下的SQL：

SELECT
	st_area ( ST_GeomFromText ( 'POLYGON ((
			112.87852 28.190856,
			112.879078 28.190562,
			112.879613 28.19024,
			112.880025 28.18993,
			112.880453 28.189525,
			112.880568 28.189433,
			112.880562 28.189324,
			112.880064 28.188912,
			112.879793 28.188709,
			112.879188 28.18857,
			112.878529 28.188726,
			112.87848 28.188768,
			112.878465 28.188882,
			112.87852 28.190856
			))', 4326 -- 指定 SRID
		) :: geography );

在数据库终端运行以下的SQL以后可以得到以下的执行结果，其结果就是面积结果。

32860.61062525073

三、结果分析

这里将统一对生成的结果进行一个简单的对比和分析。通过对比和说明，结合实际的情况，可以看到那种计算方法更加接近实际情况。

1、不同算法结果对比

在应用程序中我们使用6组数据分别调用Geotools和GeographicLib-Java库，其中的一个小学我们还使用了PostGIS的空间函数直接求解的方式。调用代码如下：

public static void main(String[] args) throws NoSuchAuthorityCodeException, FactoryException, Exception {
	//计算梅溪湖的AOI 面积
    Polygon polygon = convertAoi2Polygon(MXH_AOI);
    System.out.println("梅溪湖的球面面积为：" + calculateExactGeodeticArea(polygon));
    System.out.println("m2 " + calculateGeodeticArea(polygon, CRS.decode("EPSG:4326")));
    System.out.println("*************************************************************");
    polygon = convertAoi2Polygon(MXQX_AOI);
    System.out.println("梅溪青秀的球面面积为：" + calculateExactGeodeticArea(polygon));
    System.out.println("m2 " + calculateGeodeticArea(polygon, CRS.decode("EPSG:4326")));
    System.out.println("*************************************************************");
    polygon = convertAoi2Polygon(YMGY_AOI);
    System.out.println("杨梅公园的球面面积为：" + calculateExactGeodeticArea(polygon));
    System.out.println("m2 " + calculateGeodeticArea(polygon, CRS.decode("EPSG:4326")));
    System.out.println("**************************************************************");
    polygon = convertAoi2Polygon(BBG_AOI);
    System.out.println("梅溪新天地的球面面积为：" + calculateExactGeodeticArea(polygon));
    System.out.println("m2 " + calculateGeodeticArea(polygon, CRS.decode("EPSG:4326")));
    System.out.println("************************************************************");
    polygon = convertAoi2Polygon(THL_AOI);
    System.out.println("桃花岭风景区的球面面积为：" + calculateExactGeodeticArea(polygon));
    System.out.println("m2 " + calculateGeodeticArea(polygon, CRS.decode("EPSG:4326")));
    System.out.println("*************************************************************");
    polygon = convertAoi2Polygon(YYSX_AOI);
    System.out.println("实验小学的球面面积为：" + calculateExactGeodeticArea(polygon));
    System.out.println("m2 " + calculateGeodeticArea(polygon, CRS.decode("EPSG:4326")));
}

程序运行结果如下：

这里把面积的结果整理成以下表格。请注意，在以下的表格中，面积的单位是平方米：

序号	AOI名称	GeographicLib	Geotools
1	梅溪湖风景区	3003365.9397251983	3007843.6566917896
2	梅溪青秀	33133.088184464876	33182.49462517351
3	杨梅公园	20309.3465869762	20339.625071160495
4	步步高梅溪新天地	234694.91455752775	235044.35333895683
5	桃花岭风景区	5551323.14381226	5559616.194570124
6	岳麓区实验小学	32884.20975263743	32933.23007472232

通过以上的表格可以看到，两种不同的计算方式的结果还是有一些差异，不过这些差异并不大。

2、与互联网AOI对比

以桃花岭公园为例，桃花岭景区位于湖南长沙岳麓区，是岳麓山风景名胜区的重要组成部分。桃花岭公园规划面积约4360亩。景区内包括洪寺庵水库、腊八寺遗址等重要景点。

将4360亩换算成平方米，大约是3583333.33平方米，实际上我们算出来的面积有5551323.14381226，当然我们通过高德AOI获取的数据与原始的面积可能存在一些出入，所以导致了在计算与互联网公式面积的一个差异，当然如果有更准去的数据，也可以在评论区留言修正。

3、与天地图测面对比

在上面的小节中，可能出现实际的面积与计算的结果有较大的出入的问题，接下来我们以某小学等为例，使用天地图的测面积方法来实现面积的预估，看与实际的差距有多大，首先来看下天地图中的面积测算。该小学的面积大约是3.46公顷，合34500平方米。

来看后台的计算的结果分别是：32884.20975263743（GeographicLib）和32933.23007472232（Geotools），再来看PostGIS中的计算结果是：32860.61062525073，可以看到对于该小学的面积计算，三者都是比较接近的。同时三者与实际的面积也相差不大。

四、总结

以上就是本文的主要内容。通过实验的方式可以看到，使用PostGIS生成的面积最小，而GeographicLib其次，而使用Geotools的生成是最大的。当然，这与Geotools中的地球半径等客观因素油酸。总之，基于Java和PostGIS的AOI面数据球面面积计算是一个复杂但重要的任务。通过结合Java的编程能力和PostGIS的空间数据处理能力，可以构建一个高效、准确的球面面积计算系统。在实际应用中，需要考虑多种因素，如数据精度、投影方式、地球模型等，并通过多种方法进行验证，以确保计算结果的准确性。通过模块化设计和JDBC集成，可以提高系统的灵活性和可维护性，满足不同应用场景的需求。行文仓促，定有不足之处，欢迎各位朋友在评论区批评指正，不胜感激。