WGS84坐标系下PHP非矩形坐标区域面积计算失效求助

问题根源分析 #

你之前的两种实现都基于平面矩形面积逻辑，直接用水平/垂直距离相乘，但WGS84是椭球坐标系，经纬度围成的区域本质是球面梯形，而非平面矩形：

• 第一种方法仅用左上角点的经度差距离计算水平长度，当区域南北跨度大时，不同纬度的经度周长差异显著，单一纬度的水平距离无法代表整个区域的平均宽度；
• 第二种方法用平均纬度近似经度距离，但对于南北拉伸较大的区域，这种近似会产生明显误差。

准确解决方案（基于PHPCoord） #

利用PHPCoord内置的椭球多边形面积计算能力，直接处理经纬度围成的闭合区域，自动考虑地球曲率，适配任意长宽比：

代码实现 #

use PHPCoord\GeographicPoint;
use PHPCoord\Geographic3D;
use PHPCoord\Degree;
use PHPCoord\Metre;
use PHPCoord\Polygon;

public static function calculateAccurateArea(float $latitude1, float $longitude1, float $latitude2, float $longitude2): float
{
    // 初始化WGS84坐标系
    $wgs84CRS = Geographic3D::fromSRID(Geographic3D::EPSG_WGS_84);
    
    // 按顺时针顺序构造矩形的四个顶点
    $topLeft = GeographicPoint::create($wgs84CRS, new Degree($latitude1), new Degree($longitude1), new Metre(0));
    $topRight = GeographicPoint::create($wgs84CRS, new Degree($latitude1), new Degree($longitude2), new Metre(0));
    $bottomRight = GeographicPoint::create($wgs84CRS, new Degree($latitude2), new Degree($longitude2), new Metre(0));
    $bottomLeft = GeographicPoint::create($wgs84CRS, new Degree($latitude2), new Degree($longitude1), new Metre(0));
    
    // 创建闭合多边形
    $rectanglePolygon = Polygon::create([$topLeft, $topRight, $bottomRight, $bottomLeft]);
    
    // 计算面积并转换为平方千米（默认返回单位为平方米）
    $areaSquareMetres = $rectanglePolygon->calculateArea()->getValue();
    return $areaSquareMetres / 1000000;
}

方案优势 #

1. 原生适配椭球模型：PHPCoord的calculateArea方法基于WGS84椭球的官方参数（长半轴6378137m、扁率1/298.257223563）计算，而非简化的球体近似；
2. 无近似误差：直接计算闭合球面梯形的面积，不管区域是正方形还是极端拉伸的矩形，都能精准反映实际面积；
3. 无需手动转换：不需要额外做投影转换（如UTM），直接处理原始WGS84坐标，适配全球任意区域。

备选方案（基于proj4php投影转换） #

如果需要用平面坐标计算，可将WGS84坐标转换为对应区域的UTM平面投影，再计算矩形面积：

use proj4php\Proj4php;
use proj4php\Proj;
use proj4php\Point;

public static function calculateAreaViaUTM(float $lat1, float $lon1, float $lat2, float $lon2): float
{
    $proj4 = new Proj4php();
    // WGS84坐标系
    $wgs84 = new Proj('EPSG:4326', $proj4);
    // 自动匹配UTM投影带（示例：按左上角坐标确定带号）
    $utmZone = (int)(($lon1 + 180) / 6) + 1;
    $utmProj = new Proj("EPSG:326{$utmZone}", $proj4);
    
    // 转换四个顶点到UTM平面
    $p1 = $proj4->transform($wgs84, $utmProj, new Point($lon1, $lat1));
    $p2 = $proj4->transform($wgs84, $utmProj, new Point($lon2, $lat1));
    $p3 = $proj4->transform($wgs84, $utmProj, new Point($lon2, $lat2));
    $p4 = $proj4->transform($wgs84, $utmProj, new Point($lon1, $lat2));
    
    // 计算平面矩形面积（单位：平方米，转平方千米）
    $width = abs($p2->x - $p1->x);
    $height = abs($p4->y - $p1->y);
    return ($width * $height) / 1000000;
}

注：此方案需确保区域完全在一个UTM带内，跨带时需拆分处理。

内容的提问来源于stack exchange，提问作者Robert Blasco Villarroya