基于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：

SELECTst_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集成，可以提高系统的灵活性和可维护性，满足不同应用场景的需求。行文仓促，定有不足之处，欢迎各位朋友在评论区批评指正，不胜感激。