GO.difference bug report

[gaelforget]

In Climatology.jl we use GO.difference to calculate overlaps and differences between regions where sea surface temperature exceeds certain thresholds (see Marine Heat Waves).

However, we find that the behavior of GO.difference is not robust. In some cases it works but other times it errors out for no obvious reason. As long as that is the case, I cannot recommend adopting this method for tracking heatwaves.

Help would be greatly appreciated. Is this a well known bug? Is there hope it could be fixed quickly?

Code to reproduce error :

1. setup

using Climatology, NCDatasets, CairoMakie
import GeometryOps as GO
import GeoInterface as GI

begin
    path_OISST_stats=Climatology.downloads.OISST_stats_download()
	file_climatology=joinpath(path_OISST_stats,"OISST_mean_monthly_1992_2011.nc")

    input_path=joinpath(tempdir(),"sst_dev")
    files_list=SST_processing.download_files(path=input_path,short_demo=true)
    file_sst=joinpath(input_path,"202509","oisst-avhrr-v02r01.20250915_preliminary.nc")
    sst_anom=SST_FILES.read_map(variable="anom_recompute",file=file_sst,file_climatology=file_climatology)
end

begin
	var0=deepcopy(sst_anom); var0[ismissing.(var0)].=NaN
	
	min0=2.0; max0=Inf; min1=1.0
	pol0=GO.polygonize(0.125:0.25:359.875,-89.875:0.25:89.875, min0 .< var0 .< max0)
	pol1=GO.polygonize(0.125:0.25:359.875,-89.875:0.25:89.875, min1 .< var0 .< max0)
end

begin
    i0=9 ; i1=2;
    plot(pol1); plot!(pol0);
    p1_geom=GI.getgeom(pol1)[i1];
    p0_geom=GI.getgeom(pol0)[i0];
    plot!(p1_geom); plot!(p0_geom)
    current_figure()
end

should create the following plot

2. to reproduce error

GO.difference(p1_geom, p0_geom; target = GI.PolygonTrait())

3. error message

[rafaqz]

Seems its possible to have a zero as curr.neighbor here:

https://github.com/JuliaGeo/GeometryOps.jl/blob/main/src/methods/clipping/clipping_processor.jl#L986-L988

We need to handle this/work out how it occurs @asinghvi17 .

[rafaqz]

@gaelforget

Is this a well known bug?

Well known bugs would have issues in the queue here, there is nowhere else to know them ;)

As long as that is the case, I cannot recommend adopting this method for tracking heatwaves.
Is there hope it could be fixed quickly?

We have jobs too, this is a volunteer effort. It's best to keep this kind of pressure out of github issues unless you are financing or contributing to a project.

You could speed up the process by checking if any .neighbor values are 0 (on that line where the error occurs) and see if breaking on that case fixes the problem. We would just would have to do the same thing.

(Its also much faster if your MWE really is minimal, e.g. define a geom that breaks GO.difference without all the packages)

[gaelforget]

Seems its possible to have a zero as curr.neighbor here:

https://github.com/JuliaGeo/GeometryOps.jl/blob/main/src/methods/clipping/clipping_processor.jl#L986-L988

We need to handle this/work out how it occurs @asinghvi17 .

Not sure if this helps but I ended up digging a bit deeper, put a show statement in _trace_polynodes, and got this :

[gaelforget]

@gaelforget

Is this a well known bug?

Well known bugs would have issues in the queue here, there is nowhere else to know them ;)

Could have been that the class of methods chosen here is known to fail under known conditions. I saw the names of a couple book authors (?) dropped in the code somewhere but dont know this literature, or the theory being the algo, etc.

I did not know what to search for in the issue thread. Assumed someone might link relevant issues PRs or whatever if they existed. From your response it sounds like the underlying issue is not discussed elsewhere in the repo.

But it could also have been that the issue what discussed at length in slack 6 months ago, or in one of the many other apps that folks like to use to discuss these things.

As long as that is the case, I cannot recommend adopting this method for tracking heatwaves.
Is there hope it could be fixed quickly?

Thanks for the feedback.

First sentence is a statement of fact -- when people ask me, I just cannot recommend them using this package for that goal now. Hopefully at some point in the future you'll change my mind. Does not mean the package cannot be used for fun / exploration.

wrt the question, I was hoping for an answer like : yes, it's a matter of adding an if statement in the right place, or something like that; or, no a more general algorithm does not exist, no it would be prohibitive to implement, etc. In retrospect I probably should have just asked "when can this be fixed?" but to me that seemed less gentle.

We have jobs too, this is a volunteer effort. It's best to keep this kind of pressure out of github issues unless you are financing or contributing to a project.

I hear you. Everyone is oversubscribed, has too many deadlines, and is constantly asked to do more with less. No-one is funded exclusively/explicitely to do open source code either. We all do what we can and I definitely appreciate anyone else's efforts even though I might not always express it enough or clearly.

I have contributed to this package btw, and am not saying I wont do that again. As an academic, the question of whether I can adopt software to build something bigger often come to "how responsive are the devs? does it work now? will it work in timely fashion? is there reliable user support?"

You could speed up the process by checking if any .neighbor values are 0 (on that line where the error occurs) and see if breaking on that case fixes the problem. We would just would have to do the same thing.

No idea what .neighbor is or what it means for it to be 0 or not. Dont know who "we" is but the difference between devs and users is typically that I dont know the code internals. Makes a big difference in the amount of time and effort it'd require us to achieve the same task.

Thank you very much for the quick response btw

(Its also much faster if your MWE really is minimal, e.g. define a geom that breaks GO.difference without all the packages)

Not sure what you mean / how to do that. I took a bunch of time after work last night to make it easy to reproduce the error and write up this issue explaining what I can make out.

I thought this was a useful and nice way to help, and about as much as I had time to do.

The geom in the example breaks GO.difference and there is no other package involved.

My guess is you would have liked an example that does not involve downloading and reading a netcdf file. Seems very hard to do.

ps. sorry for all the English text but I felt compelled to explain myself

[rafaqz]

No idea what .neighbor is or what it means for it to be 0 or not.

On the line from your error that I linked, and the message is specifically saying it's indexing with zero! I don't know any more than that

[gaelforget]

No idea what .neighbor is or what it means for it to be 0 or not.

On the line from your error that I linked, and the message is specifically saying it's indexing with zero! I don't know any more than that

No worries.

I see you looped in @asinghvi17 -- thanks! Seems that @skygering and @LanaLubecke have also contributed to these functionalities, so I am looping them in here too.

[asinghvi17]

Hey! I'm currently on vacation but could look into this on one of the innumerable train rides.

I have a nice error struct that encodes all the information, seems like it doesn't catch these kinds of errors though so we'd have to extend it to that.

I think what @rafaqz meant in his request, is that it would be easier to have the polygons as a GeoJSON or something we could pick up and directly open, without the setup code. I might have the error write a temp geojson that you could copy and paste as well, at some point.

But I can probably replicate from this and will reach out to you if that proves impossible.

[gaelforget]

Hey! I'm currently on vacation but could look into this on one of the innumerable train rides.

Nice!

I have a nice error struct that encodes all the information, seems like it doesn't catch these kinds of errors though so we'd have to extend it to that.

I think what @rafaqz meant in his request, is that it would be easier to have the polygons as a GeoJSON or something we could pick up and directly open, without the setup code. I might have the error write a temp geojson that you could copy and paste as well, at some point.

Thanks for clarifying. Definitely helped.

But I can probably replicate from this and will reach out to you if that proves impossible.

From GeoJSON.write("p1_geom.json",p1_geom) (which I believe is what you meant for me to do) :

p0_geom.json
p1_geom.json

[LanaLubecke]

Hello everyone .neighbor is 0 if the node is not an intersection point between two polygons. When you are taking the difference of two polygons, the clipping algorithm generates a list of all the vertices and intersections points in polygon “A” and another list for those in polygon “B”. If the node you are considering is an intersection point between polygons A and B, and it is a node in "a_list," .neighbor gives you its index in the list of nodes for the other polygon (b_list). Please let me know if I can clarify anything! Best, Lana Lubecke Get Outlook for iOS<https://aka.ms/o0ukef>

…

________________________________ From: Gael Forget ***@***.***> Sent: Tuesday, September 23, 2025 6:18:14 AM To: JuliaGeo/GeometryOps.jl ***@***.***> Cc: Lubecke, Lana C. (Lana) ***@***.***>; Mention ***@***.***> Subject: Re: [JuliaGeo/GeometryOps.jl] GO.difference bug report (Issue #335) [https://avatars.githubusercontent.com/u/20276764?s=20&v=4]gaelforget left a comment (JuliaGeo/GeometryOps.jl#335)<#335 (comment)> No idea what .neighbor is or what it means for it to be 0 or not. On the line from your error that I linked, and the message is specifically saying it's indexing with zero! I don't know any more than that No worries. I see you looped in @asinghvi17<https://github.com/asinghvi17> -- thanks! Seems that @skygering<https://github.com/skygering> and @LanaLubecke<https://github.com/LanaLubecke> have also contributed to these functionalities, so I am looping them in here too. — Reply to this email directly, view it on GitHub<#335 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BAWG35ERH4P3A475DHJCXHT3UFCBNAVCNFSM6AAAAACHGOIAXOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTGMRTHE4DIOJVGA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[asinghvi17]

One problem is that polygonize is generating invalid polygons.

Here is setup code:

using GeoJSON
import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG

p0 = GeoJSON.read("p0_geom.json")
p1 = GeoJSON.read("p1_geom.json")

julia> LG.isValidReason(p0 |> GI.convert(LG))
"Ring Self-intersection[78.25 72.25]"

julia> LG.isValidReason(p1 |> GI.convert(LG))
"Ring Self-intersection[78.25 72.25]"

Here's a picture of raw p0. Blue is exterior, red is holes.

You can 'correct' this with LG:

julia> GI.nhole(p0)
12

julia> LG.makeValid(p0 |> GI.convert(LG)) |> GI.nhole
16

Here's what the fixed one looks like, again blue is exterior and red is holes. Now the self intersection is forming a hole.

Once this is fixed, then GeometryOps does fine here.

So for @gaelforget's use case if we fix polygonize then all is good. But we should definitely either detect this or throw a better error here, ideally both.

[asinghvi17]

Just to clarify this is also a bug in difference that it doesn't error when you feed it an invalid polygon, but the usecase can probably be enabled by just generating correct polygons instead.

[asinghvi17]

Also, you can manually enable tree acceleration to get an ~16x speedup on this problem:

import LibGEOS as LG, GeometryOps as GO
p0f = LG.makeValid(GI.convert(LG, pol0)) |> GO.tuples
p1f = LG.makeValid(GI.convert(LG, pol1)) |> GO.tuples


tree_alg = GO.FosterHormannClipping(
    GO.SingleNaturalTree()
)
res1 = GO.difference(p0f, p1f; target = GI.PolygonTrait())
res2 = GO.difference(tree_alg, p0f, p1f; target = GI.PolygonTrait())
@assert isempty(res1) == isempty(res2) || 
        GO.equals(only(res1), only(res2))

using Chairmarks
# naive O(nm) - 25 ms
@be GO.difference($p0f, $p1f; target = $(GI.PolygonTrait()))
# smart tree alg - 1.5 ms
@be GO.difference($tree_alg, $p0f, $p1f; target = $(GI.PolygonTrait()))

If you reverse the order of the input arguments, then you only get a 3x speedup because it seems to take 7ms to actually assemble the polygons. But it's still quite significant.