GIS中图斑椭球面积的计算

我们都知道投影就是将球面通过一定的规则展开成平面，方便浏览和计算，投影就会产生变化，无论是角度、面积还是长度。在实际工作中，经常要计算图斑面积，这里简单总结一下常用的计算方法。

面积量算

先来看一下ArcMap中不同坐标系统通过面积量算得到的结果。

情况一：无空间参考，在空白地图中加载一个没有定义坐标系统的图斑（虽然没有坐标系统，但有坐标值，只不过坐标单位是未知单位），在加入时会弹出"未知的空间参考"的警告。

使用测量工具能够量算距离和面积，虽然量测值，但是量测单位为未知，当然也就不能切换距离或面积单位了。

情况二：地理坐标。我们将上述图斑文件定义正确的地理坐标，重新打开一个空白地图添加定义空间参考后的图斑文件，发现可以测量距离，坐标单位可以由度（数据本身的单位）切换成平面单位，但无法量测面积。

情况三：投影坐标。我们将上述定义了地理坐标的文件投影成合适的投影坐标（38带），就可以正常量取距离和面积了。

情况四：不同投影坐标。我们再将上述定义了地理坐标的文件投影另外一个投影带（30带），对比可看出它们二者的形变。

通过计算几何，可看到它们的面积差距是有差距的。

计算几何和通过!shape.Area!计算的结果都一样，都是投影面积。

一般地，采用高斯3度投影，相邻的两个投影带计算的面积差距可能几十到几百平方不等。因此，这个面积只能是一个粗略值，在一些对面积精度要求比较高的项目中（如承包地）则不能采用投影面积。

椭球面积

那么，更准确一点的就是椭球面积，这里不讲椭球面积的计算公式，在二调项目时，国土资源部下发了《关于统一图幅理论面积与图斑椭球面积计算要求》的通知，对图幅理论面积与图斑椭球面积计算公式进行了细化，明确了面积计算方法，统一了公式中的有关参数，有兴趣的可以自行查询和编程实现。

在ArcMap使用!shape.geodesicArea!即可计算，可以看出，面积1（38带）更接近椭球面积，而面积2（30带）则差距很大，因此选择一个正确的投影非常重要。

值得注意的：

（1）同一个数据，在其地理坐标系统和投影坐标系统下计算的椭球面积有细微差距，原因是投影坐标计算椭球面积时，要先将投影转换成地理坐标模型，这个转换会有误差。

（2）如果需要生成指定单位的椭球面积，可以在后面@单位，如!shape.geodesicArea@SQUAREMETERS!，表示的单位是平方米，如果不定，则是坐标系统的默认单位。

（3）未定义任何坐标系统的数据无法计算椭球面积。

代码实现

使用Python代码，利用ArcPy计算：

# 投影面积计算
arcpy.CalculateField_management(fc,"ProjectArea","!shape.Area!","PYTHON_9.3")
# 椭球面积计算
arcpy.CalculateField_management(fc,"GeoArea","!shape.geodesicArea@SQUAREKILOMETERS!","PYTHON_9.3")

使用C#代码，利用ArcObject接口计算：

线接口：IPolycurveGeodetic.LengthGeodetic

面接口：IAreaGeodetic.AreaGeodetic

private void GetArea(IFeatureClass featureClass)
{
    IFeatureCursor featureCursor = featureClass.Search(null, true);
    IFeature feature = featureCursor.NextFeature();
    Console.WriteLine($@"序号          投影面积         正截面          测地线        等角航线        大椭圆面积");
    while (feature != null)
    {
        IArea area = (IArea)feature.Shape;
        IAreaGeodetic areaGeodetic = (IAreaGeodetic)feature.Shape;
        ILinearUnit linearUnit = new LinearUnitClass();
        double geodeticArea1 = areaGeodetic.AreaGeodetic[esriGeodeticType.esriGeodeticTypeNormalSection, linearUnit];
        double geodeticArea2 = areaGeodetic.AreaGeodetic[esriGeodeticType.esriGeodeticTypeGeodesic, linearUnit];
        double geodeticArea3 = areaGeodetic.AreaGeodetic[esriGeodeticType.esriGeodeticTypeLoxodrome, linearUnit];
        double geodeticArea4 = areaGeodetic.AreaGeodetic[esriGeodeticType.esriGeodeticTypeGreatElliptic, linearUnit];
        Console.WriteLine($@"{feature.OID}           {area.Area}         {geodeticArea1}         {geodeticArea2}         {geodeticArea3}           {geodeticArea4}");
        feature = featureCursor.NextFeature();
    }
}