Could you please help me understand the issue with H3 geospatial indexing?
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)
I'm expecting to get indices like these ones, which are located inside of the geojson polygon above:
But instead, h3.polyfill
returns indices that are “flipped” by 90 degrees like these:
All works fine for other geojsons that don't span Antarctica.. I’m using Python 3.10.7 and H3 3.7.4.
I would appreciate any hints.
Upd.
I used geo_json_conformant=True
parameter and it flipped indices back. But it seems not all resolution 3 indices were generated and my expected indices are not in the list. On the image generated indices are in blue and expected are in red.
Upd 2
Following the suggestion from @nrabinowitz, I triangulated the original polygon from Pole and then polyfilled resulting "slices". Works perfectly fine, all missing indices are in place. result
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))
H3 uses a Cartesian model for polygons used in polyfill
, not a spherical model, so poles are slightly challenging. The other issue is that we assume the smaller polygon when some of the arcs are greater than 180 degrees, which is frequently the case very near the poles.
See the suggested workaround in this related issue:
I think the best workaround here is to slice up polygons that contain a pole. I think the simplest version of this, which ought to work, is to make triangles, one for each pair of vertexes, with the third vertex being the pole itself.
The other participant in that discussion glossed this as "slicing the polygon up like a pizza," which I thought was very descriptive.
There's a demo of this approach in this Observable notebook.