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.

[juliohm]

I will think about this specific method, to see if it makes sense to add a method to sideof with MultiRing or a method to in with Polygon with holes. I thought we already had one.

[juliohm]

The problem is that we are restricting the method currently to Euclidean polygons:

Meshes.jl/src/predicates/in.jl

Line 113 in e959de5

function Base.in(point::Point, poly::Polygon{𝔼{2}})

Should this method get relaxed to also support curved polygons over the ellipsoid? I understand that your use case involves points with LatLon coordinates.

[juliohm]

We need to think carefully about these geodesic variants later. Maybe the easiest thing to do now is to Proj the dataset, run the predicate in a flat space, and then undo the Proj if necessary. What do you think?

[juliohm]

This is the current situation:

point in dom already exists as you said:

Meshes.jl/src/predicates/in.jl

Line 160 in 097754c

Base.in(p::Point, d::Domain) = any(e -> p ∈ e, d)

point in poly is only implemented for Euclidean polygons:

Meshes.jl/src/predicates/in.jl

Line 113 in e959de5

function Base.in(point::Point, poly::Polygon{𝔼{2}})

When you call point in dom with point and dom over the ellipsoid, the fallback of dom calls the fallback of Geometry:

Meshes.jl/src/predicates/in.jl

Line 10 in 097754c

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

We need to figure out the correct solution here. We plan to look into geodesic geometry after we finish polishing the Euclidean geometry algorithms.

[disberd]

I see, on my side it is fine to just keep getting the error for Point(LatLon) until we have some geodesic computations in place to be consistent.

A potential fallback till then is to call flat on the LatLon coords to make them Cartesian2D (which is somehow alredy what happens when just do p in poly at the moment if poly does not have holes as far I understand.

In any case, I understand if you want this to error for LatLon (though it should probably then error for every type of polygon (even without holes) for consistency.

(For my use case I can just flatten the coords myself before checking for inclusion, knowing I am doing approximation since distance on E{2} is not the same as distance on ellipsoid)

[juliohm]

An alternative method consists of converting LatLon to PlateCarree, which is basically the naive approach that treats 1 deg = 1 m. It is what happens when people plot LatLon coordinates in a 2D Cartesian system:

flatpoint = point |> Proj(PlateCarree)
flattable = geotable |> Proj(PlateCarree)

flatpoint in flattable.geometry

In any case, I understand if you want this to error for LatLon (though it should probably then error for every type of polygon (even without holes) for consistency.

That is a good point. We should probably throw an error when a geodesic geometry is given as input to the predicate. And then, when we are ready to start the work on geodesic algorithms, we can add methods one by one.