Algorithms from geodesic geometry

[disberd]

I am trying to do some benchmarks to improve the performance of sideof (see #988), and I found method errors very early with the following example, probably caused by missing methods for sideof:

using Meshes, GeoArtifacts

dmn = NaturalEarth.get("admin_0_countries", 110).geometry
point_nemo = Point(LatLon(−48.8767,−123.3933))

point_nemo in dmn

This throws a first error

ERROR: TypeError: in typeassert, expected Integer, got a value of type Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}
Stacktrace:
 [1] _similar_shape(itr::Base.Generator{Point{🌐, GeodeticLatLon{…}}, Meshes.var"#309#310"{Box{…}}}, ::Base.HasLength)
   @ Base ./array.jl:652
 [2] collect(itr::Base.Generator{Point{🌐, GeodeticLatLon{…}}, Meshes.var"#309#310"{Box{…}}})
   @ Base ./array.jl:779
 [3] sideof(points::Point{🌐, GeodeticLatLon{…}}, object::MultiRing{🌐, GeodeticLatLon{…}, Ring{…}})
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/sideof.jl:163
 [4] in
   @ ~/.julia/packages/Meshes/FJ9DT/src/predicates/in.jl:10 [inlined]
 [5] (::Meshes.var"#280#281"{Point{…}})(e::PolyArea{🌐, GeodeticLatLon{…}, Ring{…}, Vector{…}})
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/predicates/in.jl:160
 [6] _any(f::Meshes.var"#280#281"{Point{…}}, itr::GeometrySet{🌐, GeodeticLatLon{…}, Geometry{…}}, ::Colon)
   @ Base ./reduce.jl:1237
 [7] any
   @ ./reduce.jl:1228 [inlined]
 [8] in(p::Point{🌐, GeodeticLatLon{…}}, d::GeometrySet{🌐, GeodeticLatLon{…}, Geometry{…}})
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/predicates/in.jl:160
 [9] top-level scope
   @ REPL[26]:1
Some type information was truncated. Use `show(err)` to see complete types.

which is caused by this fallback method not working when a single Point is provided:

Meshes.jl/src/sideof.jl

Lines 161 to 168 in c0c9c2d

	function sideof(points, object::GeometryOrDomain)
	bbox = boundingbox(object)
	isin = [point ∈ bbox for point in points]
	inds = findall(isin)
	side = fill(OUT, length(isin))
	side[inds] .= sidewithinbox(collectat(points, inds), object)
	side
	end

This can be easily fixed by creating a method for the single point that just wraps it in an iterable (here in the REPL for simplicity)

julia> @eval Meshes sideof(p::Point, obj::GeometryOrDomain) = sideof((p,), obj)

After fixing that method, another error appears though

 julia> point_nemo in dmn                                                                                            
ERROR: MethodError: no method matching sidewithinbox(::Vector{Point{…}}, ::MultiRing{🌐, GeodeticLatLon{…}, Ring{…}})                                                                                                                   The function `sidewithinbox` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  sidewithinbox(::Any, ::Mesh)
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/sideof.jl:172
  sidewithinbox(::Any, ::Ring)
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/sideof.jl:170

Stacktrace:
 [1] sideof(points::Tuple{Point{…}}, object::MultiRing{🌐, GeodeticLatLon{…}, Ring{…}})
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/sideof.jl:166
 [2] sideof(p::Point{🌐, GeodeticLatLon{…}}, obj::MultiRing{🌐, GeodeticLatLon{…}, Ring{…}})
   @ Meshes ./REPL[27]:1
 [3] in
   @ ~/.julia/packages/Meshes/FJ9DT/src/predicates/in.jl:10 [inlined]
 [4] (::Meshes.var"#280#281"{Point{…}})(e::PolyArea{🌐, GeodeticLatLon{…}, Ring{…}, Vector{…}})
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/predicates/in.jl:160
 [5] _any(f::Meshes.var"#280#281"{Point{…}}, itr::GeometrySet{🌐, GeodeticLatLon{…}, Geometry{…}}, ::Colon)
   @ Base ./reduce.jl:1237
 [6] any
   @ ./reduce.jl:1228 [inlined]
 [7] in(p::Point{🌐, GeodeticLatLon{…}}, d::GeometrySet{🌐, GeodeticLatLon{…}, Geometry{…}})
   @ Meshes ~/.julia/packages/Meshes/FJ9DT/src/predicates/in.jl:160
 [8] top-level scope
   @ REPL[28]:1
Some type information was truncated. Use `show(err)` to see complete types.

which again seems to be caused by some missing method implemented

[juliohm]

Thanks for catching these missing methods. Please feel free to submit a PR. We can merge and release a patch quickly.

[juliohm]

@disberd please keep in mind that some methods aren't supposed to exist, for example point_nemo in dmn only makes sense when dmn is a surface mesh. In this case what you really want is broadcast: point_nemo .in dmn

[juliohm]

Are there any other methods missing after this broadcast is applied?

[disberd]

@juliohm my use case is to verify whether a given point is inside any of a subset of countries (so a subdomain out of the full dmn). Given that in(::Point, ::Domain) dispatches to any(e -> p in e, d) I believe that to be exactly what I want to achieve. I am not interested in a vector of results of in with every country (any is also potentially faster as it exits early instead of computing in for every polygon if one of the early polygons returns true).

Is there a better approach to achieve what I want above?

That being said point_nemo .in dmn synthax is not parsed correctly, while point_nemo .∈ dmn just gives the same errors above

[juliohm]

my use case is to verify whether a given point is inside any of a subset of countries (so a subdomain out of the full dmn).

any(geom -> point ∈ geom, dom)

should do it right? The only methods for sideof are (point, line), (point, ring) and (point, mesh) with a surface mesh because these objects split the space into regions (e.g., left and right half-spaces, interior and exterior volume)

[juliohm]

That being said point_nemo .in dmn synthax is not parsed correctly,

It works for me here:

triangles = rand(Triangle, 100) |> Shadow("xy")

point = rand(Point) |> Shadow("xy")

point .∈ triangles

[disberd]

my use case is to verify whether a given point is inside any of a subset of countries (so a subdomain out of the full dmn).

any(geom -> point ∈ geom, dom)

should do it right? The only methods for sideof are (point, line), (point, ring) and (point, mesh) with a surface mesh because these objects split the space into regions (e.g., left and right half-spaces, interior and exterior volume)

But this is already what happens by default with point in dom, so isn't it just easier to do point in dom? (I was not calling sideof directly above)

[disberd]

That being said point_nemo .in dmn synthax is not parsed correctly,

It works for me here:

triangles = rand(Triangle, 100) |> Shadow("xy")

point = rand(Point) |> Shadow("xy")

point .∈ triangles

Yes the .∈ works, while .in is not parsed by julia, the problem of the errors above just come from PolyAreas with both an inner and outer ring (so I guess polygons with holes). Those are internally transformed into this call:

Meshes.jl/src/predicates/in.jl

Line 10 in e959de5

Base.in(p::Point, g::Geometry) = sideof(p, boundary(g)) == IN

which apparently transforms the PolyArea into a MultiRing (inside boundary) and creates the issue of sideof

[juliohm]

Oh I see, so we are missing a method of sideof with MultiRing. Also this generic fallback should probably do != OUT instead, given that ON is also considered and geometries are closed.