Investigate failing rasterization

[moradology]

This combination of polygon and raster extent fails during rasterization with default rasterizer settings despite what appears to be a valid geometry. It would be nice to get better error messages/work around cases like this if possible

Polygon

POLYGON ((-156.7523879306299 21.318054818524036, -156.75149525021493 21.320234352774555, -156.75061988864124 21.322420900278683, -156.74976190100162 21.32461432342134, -156.74892134129536 21.326814484154717, -156.7480982624249 21.32902124400696, -156.74729271619248 21.331234464090887, -156.74650475329688 21.333454005112728, -156.74573442333022 21.335679727380885, -156.7449817747749 21.337911490814736, -156.7442468550005 21.340149154953448, -156.74352971026076 21.342392578964805, -156.74283038569072 21.344641621654088, -156.7421489253039 21.34689614147295, -156.74148537198943 21.349155996528335, -156.7408397675095 21.351421044591397, -156.7402121524966 21.35369114310646, -156.73960256645108 21.355966149199986, -156.7390110477385 21.358245919689562, -156.73843763358738 21.36053031109293, -156.7378823600868 21.362819179636993, -156.73734526218402 21.365112381266886, -156.73682637368253 21.367409771655023, -156.73632572723966 21.36971120621019, -156.73584335436468 21.372016540086655, -156.73537928541677 21.374325628193258, -156.7349335496031 21.376638325202574, -156.734506174977 21.378954485560033, -156.73409718843627 21.381273963493094, -156.73370661572127 21.383596613020423, -156.7333344814136 21.385922287961066, -156.73298080893434 21.388250841943666, -156.73264562054257 21.390582128415662, -156.73232893733413 21.39291600065252, -156.73203077924018 21.395252311766967, -156.7317511650259 21.39759091471823, -156.73149011228946 21.399931662321297, -156.7312476374607 21.40227440725618, -156.73102375580038 21.40461900207718, -156.73081848139896 21.40696529922217, -156.73063182717576 21.409313151021887, -156.73046380487833 21.41166240970923, -156.73031442508147 21.414012927428537, -156.7301836971868 21.416364556244915, -156.73007162942184 21.418717148153544, -156.7299782288399 21.421070555088985, -156.72990350131937 21.423424628934505, -156.7298474515633 21.425779221531393, -156.72981008309938 21.4281341846883, -156.72979139827942 21.430489370190543, -156.72979139827942 21.432844629809455, -156.72981008309938 21.4351998153117, -156.7298474515633 21.43755477846861, -156.72990350131937 21.439909371065497, -156.7299782288399 21.442263444911017, -156.73007162942184 21.444616851846455, -156.7301836971868 21.446969443755084, -156.73031442508147 21.449321072571465, -156.73046380487833 21.451671590290772, -156.73063182717576 21.454020848978114, -156.73081848139896 21.456368700777833, -156.73102375580038 21.458714997922822, -156.7312476374607 21.461059592743823, -156.73149011228946 21.4634023376787, -156.7317511650259 21.465743085281773, -156.73203077924018 21.468081688233035, -156.73232893733413 21.470417999347482, -156.73264562054257 21.47275187158434, -156.73298080893434 21.475083158056336, -156.7333344814136 21.477411712038936, -156.73370661572127 21.47973738697958, -156.73409718843627 21.48206003650691, -156.734506174977 21.48437951443997, -156.7349335496031 21.486695674797428, -156.73537928541677 21.489008371806744, -156.73584335436468 21.491317459913347, -156.73632572723966 21.49362279378981, -156.73682637368253 21.49592422834498, -156.73734526218402 21.498221618733115, -156.7378823600868 21.500514820363005, -156.73843763358738 21.502803688907072, -156.7390110477385 21.50508808031044, -156.73960256645108 21.507367850800016, -156.7402121524966 21.509642856893542, -156.7408397675095 21.511912955408604, -156.74148537198943 21.514178003471663, -156.7421489253039 21.51643785852705, -156.74283038569072 21.518692378345914, -156.74352971026076 21.520941421035197, -156.7442468550005 21.523184845046554, -156.7449817747749 21.525422509185265, -156.74573442333022 21.527654272619117, -156.74650475329688 21.529879994887274, -156.74729271619248 21.532099535909115, -156.7480982624249 21.534312755993042, -156.74892134129536 21.536519515845285, -156.74976190100162 21.538719676578662, -156.75061988864124 21.54091309972132, -156.75149525021493 21.543099647225443, -156.7523879306299 21.545279181475966, -156.75200324049055 21.545438525349006, -156.75110930806937 21.5432559342482, -156.7502327187796 21.541066320057535, -156.74937352779125 21.538869820585084, -156.7485317891794 21.536666574072278, -156.7477075559206 21.53445671918518, -156.7469008798898 21.53224039500576, -156.7461118118569 21.530017741023162, -156.74534040148356 21.527788897124896, -156.7445866973202 21.52555400358805, -156.74385074680276 21.52331320107046, -156.74313259624992 21.521066630601865, -156.74243229086002 21.51881443357501, -156.74174987470826 21.516556751736765, -156.741085390744 21.5142937271792, -156.74043888078793 21.512025502330637, -156.73981038552958 21.509752219946687, -156.73919994452464 21.50747402310128, -156.73860759619257 21.505191055177637, -156.73803337781402 21.502903459859255, -156.73747732552872 21.50061138112088, -156.73693947433296 21.498314963219418, -156.7364198580776 21.496014350684874, -156.73591850946582 21.493709688311256, -156.73543546005104 21.491401121147447, -156.73497074023504 21.489088794488097, -156.73452437926596 21.486772853864466, -156.73409640523647 21.484453445035257, -156.73368684508205 21.482130713977465, -156.73329572457925 21.47980480687717, -156.73292306834406 21.477475870120337, -156.73256889983045 21.475144050283625, -156.73223324132877 21.472809494125134, -156.73191611396436 21.47047234857518, -156.73161753769637 21.468132760727052, -156.73133753131629 21.46579087782775, -156.73107611244694 21.463446847268717, -156.73083329754132 21.46110081657656, -156.73060910188147 21.45875293340378, -156.73040353957762 21.45640334551946, -156.73021662356732 21.45405220079996, -156.7300483656145 21.45169964721965, -156.72989877630883 21.44934583284155, -156.729767865065 21.446990905808036, -156.7296556401223 21.444635014331514, -156.72956210854377 21.442278306685086, -156.72948727621602 21.43992093119322, -156.72943114784883 21.43756303622242, -156.7293937269747 21.435204770171886, -156.72937501594882 21.43284628146417, -156.72937501594882 21.430487718535833, -156.7293937269747 21.428129229828116, -156.72943114784883 21.42577096377758, -156.72948727621602 21.42341306880678, -156.72956210854377 21.421055693314916, -156.7296556401223 21.418698985668488, -156.72976786506504 21.416343094191962, -156.72989877630883 21.41398816715845, -156.7300483656145 21.41163435278035, -156.73021662356732 21.40928179920004, -156.73040353957762 21.406930654480544, -156.73060910188147 21.40458106659622, -156.73083329754132 21.40223318342344, -156.73107611244694 21.399887152731285, -156.73133753131629 21.39754312217225, -156.73161753769637 21.39520123927295, -156.73191611396436 21.392861651424823, -156.73223324132877 21.39052450587487, -156.73256889983045 21.388189949716377, -156.73292306834406 21.385858129879665, -156.73329572457925 21.383529193122833, -156.73368684508205 21.381203286022537, -156.73409640523647 21.378880554964745, -156.73452437926596 21.376561146135536, -156.73497074023504 21.374245205511905, -156.73543546005104 21.371932878852554, -156.73591850946582 21.369624311688746, -156.7364198580776 21.367319649315128, -156.73693947433296 21.365019036780584, -156.73747732552872 21.362722618879122, -156.73803337781402 21.360430540140744, -156.73860759619257 21.358142944822365, -156.73919994452464 21.35585997689872, -156.73981038552958 21.353581780053315, -156.74043888078793 21.351308497669365, -156.741085390744 21.3490402728208, -156.74174987470826 21.346777248263237, -156.74243229086002 21.34451956642499, -156.74313259624992 21.342267369398137, -156.74385074680276 21.34002079892954, -156.7445866973202 21.337779996411953, -156.74534040148356 21.335545102875106, -156.7461118118569 21.33331625897684, -156.7469008798898 21.331093604994237, -156.7477075559206 21.328877280814822, -156.7485317891794 21.326667425927724, -156.74937352779125 21.324464179414917, -156.7502327187796 21.322267679942467, -156.75110930806937 21.320078065751797, -156.75200324049055 21.317895474650996, -156.7523879306299 21.318054818524036))

RasterExtent

GridExtent(Extent(-157.026672, 21.304191901569745, -156.69356613551687, 21.559142098430257),CellSize(3.3310586448311594E-4,6.37375492151282E-4),1000x400)

[atararaksin]

I also came across this some time ago. The trace is as follows (for @moradology's case):

Exception in thread "main" java.lang.IllegalStateException: Rasterizer encountered an error: an odd number of X axis intersections encountered via a horizontal line across the raster at y = 21.
	at geotrellis.raster.rasterize.polygon.PolygonRasterizer$.foreachCellByPolygon(PolygonRasterizer.scala:432)
	at geotrellis.raster.rasterize.Rasterizer$.foreachCellByGeometry(Rasterizer.scala:117)
	at geotrellis.raster.rasterize.Rasterizer$.foreachCellByGeometry(Rasterizer.scala:93)
	at geotrellis.raster.rasterize.Rasterizer$.rasterizeWithValue(Rasterizer.scala:67)

Looks like this indeed is a bug, because as @moradology mentioned, the polygon is valid. I've made a little investigation and, in short, I think it's a Double-related thing. Below is a more extended description of what causes the problem. I realize it looks messy, still maybe it will help; never mind if it's impossible to follow along:)

I assume that a variation of the scan-line algorithm is used, as described here: https://www.cs.rit.edu/~icss571/filling/how_to.html
There's this corner case when the scan-line crosses a vertex of the polygon:

This code seems to handle this case:

geotrellis/raster/src/main/scala/geotrellis/raster/rasterize/polygon/PolygonRasterizer.scala

Lines 106 to 107 in b01b6e2

	if (y == y1) (x1,1) // Top endpoint
	else if (y == y2) (x2,-1) // Bottom endpoint

But in our problem case, the scan-line (y=21.5) doesn't cross the vertex exactly (y=21.499999999998117). So the above-mentioned code doesn't treat the intersection point as a vertex. However, due to the imperfect precision of Double, the calculations of intersection x-coords for both edges yield the same Double. Which effectively means that the scan-line does intersect the edges in a vertex (which it doesn't). This leads to the parity sum being calculated as 0 + 0 = 0 later, while given a perfect precision it should have been either 2 or -2. As a result the point is added to rowRuns and later crashes on

geotrellis/raster/src/main/scala/geotrellis/raster/rasterize/polygon/PolygonRasterizer.scala

Line 426 in b01b6e2

if(rowRuns.length % 2 != 0) {

The most straight-forward fix could be to use BigDecimal or spire types in lineAxisIntersection, but I wonder if it's acceptable performance-wise. Another option could be to tweak these conditions:

geotrellis/raster/src/main/scala/geotrellis/raster/rasterize/polygon/PolygonRasterizer.scala

Lines 106 to 107 in b01b6e2

	if (y == y1) (x1,1) // Top endpoint
	else if (y == y2) (x2,-1) // Bottom endpoint

[pomadchin]

wow @atararaksin thank you so much for looking into it. 💯 added a bug label

If you have any ideas with fixing it with spire types / BigDecimal I'll be happy to assist you btw.

For sure as you said, there can be some performance costs so in this case these changes would have to be benchmarked :/

[atararaksin]

@pomadchin I've set up a little benchmark and yeah, using BigDecimal or spire does affect the performance (no surprise). So this will not work, but I have some other ideas, just needs a bit of a thought. I hope to come up with a PR shortly.

[atararaksin]

So basically I chose to use Rational, but to do it only for the problem points - this way the general-case performance is not affected, because for most polygons and extents this 'extra precision' is triggered zero to very few times. For a moderately large geometry set that I used for benchmarking, it was not triggered at all.

This approach may look a bit ad-hoc, but it preserves the original algorithm logic and structure, doesn't affect performance, solves the issue and is unlikely to break anything.