利用FME计算椭球面积

前言

“三调”以来，自然资源相关业务中很多地方都会用到椭球面积。学习了《2020FME博客大赛——FME在GIS中的椭球面积计算和高斯反算》之后，自己也用FME做了一个。

原理

《第三次全国国土调查技术规程》（TD/T 1055-2019）中的附录D有详细描述，截图如下：

思路

上面的公式，可以利用 PythonCaller 调用python实现；
计算需要大地坐标，可以利用 Reprojector 将投影坐标系等其他坐标系转为地理坐标系；
对于内部有空洞的图斑，可以用 DonutHoleExtractor 分离环；
StatisticsCalculator 能分组汇总面积；
FeatureJoiner能将计算得到的面积关联过来。

具体实现

套用技术规程中的公式和常数，利用Python计算椭球面积，关键代码如上图所示。

自定义转换器使用

此部分针对看完以上内容仍不会利用FME计算椭球面积的同学，高手请跳过。
这里有一个做好的转换器： 链接
由2020版FME编写，需在2020及以上版本中使用。
因需要用 Reprojector 将坐标系转为 2000国家大地坐标系，需确认原始数据有正确的坐标系才能正确计算椭球面积。

下载好后右键点击安装，回到FME Workbench中即可调用此自定义转换器

遇到的bug

FeatureJoiner这一步遇到过报错 INVALID_GEOMETRY_COORDINATE_SYSTEM 导致丢失图斑。
更换为FeatureMerger转换器后正常。

同样的数据，转换前使用CoordinateSystemSetter设置坐标系后，FeatureJoiner这一步也正常。
百思不得其解…

更新（线段长度大于70m内插点计算椭球面积）

有小伙伴反馈，技术规程提到在椭球面积计算过程中,线段长度大于70 m需要进行插值计算，计算过程怎么没有这个？
安排！

思路

Chopper能把面变成线段。
计算线段长度，除以70，向上取整，判断是否需要内插点。
利用snip()和appendLine()截取并拼接线段，完成内插点。
回到之前的思路计算椭球面积。

具体如下图：

其中完成内插点工作的PythonCaller代码块：

def input(self,feature):
    line = feature.getGeometry()
    length = line.getLength(False)
    numLines = math.ceil(length/70)
    if numLines > 1:
        len = length / numLines
        line.snip(fmeobjects.SNIP_DISTANCE, False, 0, len)
        for n, (beg, end) in enumerate([(i*len, (i+1)*len) for i in range(1,numLines)], 2):
            line2 = feature.getGeometry()
            if n == numLines:
                end = -1
            line2.snip(fmeobjects.SNIP_DISTANCE, False, beg, end)
            line.appendLine(line2)
        feature.setGeometry(line)
    self.pyoutput(feature)

之前的算法，可以计算地理坐标的图形的椭球面积。
更新后的算法，因为计算线段长度的时候没管坐标系统，只适用于投影坐标的图形。