H3 polyfill处理南极洲GeoJSON时索引异常问题求助

问题描述 #

使用H3的polyfill方法处理南极洲区域的GeoJSON多边形时出现异常，代码如下：

import h3
geo_antarctic = {"type":"Polygon","coordinates":[[[-170.63764683701507,-85.05113000000047],[-170.63764683701507,-77.89462449499929],[-63.82520589349025,-66.39564184408599],[-49.69216225292467,-77.30460454007437],[-35.16653406678777,-77.89462449499929],[-9.255954059083527,-70.29658532122083],[40.994867774038596,-68.50197979740217],[89.56411960844528,-64.94027568501143],[163.48124599227498,-67.77106116580279],[172.90327508598565,-72.42721956336818],[165.83675326570284,-77.7288586062699],[178.18462781512582,-77.47601087207454],[178.57721236069702,-85.0171471646522],[-178.63764683701507,-85.05113000000047]]]}

idx = h3.polyfill(geo_antarctic, 3)

预期得到的索引：

• 83ef9efffffffff
• 83eea4fffffffff
• 83f125fffffffff
• 83f2a4fffffffff

实际返回的索引（出现“翻转90度”的异常）：

• 836682fffffffff
• 830e59fffffffff
• 83b294fffffffff
• 836733fffffffff
• 838f0bfffffffff
• 830372fffffffff

环境：Python 3.10.7、H3 3.7.4，处理其他非南极洲区域的GeoJSON时一切正常。

第一次尝试：使用geo_json_conformant=True参数 #

添加geo_json_conformant=True参数后，索引的“翻转”问题得到解决，但仍存在遗漏：部分预期的3级索引未被生成。可视化显示生成的索引为蓝色，预期索引为红色，目标索引不在结果列表中。

最终解决方法：极点三角剖分 #

将原始多边形从南极点进行三角剖分，对每个生成的“扇形切片”单独执行polyfill，最终完美覆盖所有目标索引。实现代码如下：

import h3
import geojson

geoj = {"type":"Polygon","coordinates":[[[-178.34111242523068,-85.0207089708011],[-178.69267492523034,-77.91567194747755],[-162.52079992523068,-78.4905544838336],[-140.02079992523,-73.8248242864237],[-126.66142492523065,-73.12494935304983],[-103.10673742523004,-74.59011176731619],[-103.45829992523063,-71.07406105104535],[-83.06767492523001,-73.52840349816283],[-61.97392492523001,-64.32087770911836],[-57.052049925230655,-62.43108077917767],[-59.86454992522999,-74.77584672076205],[-39.12236242523063,-77.8418507294947],[-12.052049925230301,-70.61261893331015],[35.05732507477002,-68.52824009786191],[53.33857507476973,-65.51296841598038],[76.54170007476968,-68.39918525054024],[93.06513757477003,-64.77413134111099],[143.69013757477003,-66.08937000336596],[173.22138757477006,-70.72898413027124],[167.94795007477003,-76.26869800825351],[177.79170007476975,-77.23507678015689],[178.60169170931843,-84.94715491814792],[-178.34111242523068,-85.0207089708011]]]}

polygon_coords = geoj["coordinates"][0]
pole_coord = (0.0, -89.999)

all_indexes = set()

for i in range(len(polygon_coords)-1):
    polygon = geojson.Polygon([[pole_coord, tuple(polygon_coords[i]), tuple(polygon_coords[i+1])]])
    idxes = h3.polyfill(dict(polygon), 5, geo_json_conformant=True)
    all_indexes.update(idxes)

with open(f"./absent_polr.csv", "w") as out:
    out.write("h3_idx\r\n")
    out.write("\r\n".join(all_indexes))

内容的提问来源于stack exchange，提问作者Arty