Yes, we only have a method that checks each point at a time. There are other algorithms as others pointed out for multiple points at once. Pull requests are welcome.
I ended up making a Julia function from wiki’s algorithm pseudocode… works fast enough and no memory allocations, but its too off the ideology of Meshes.jl. would be not simple to integrate. I have a case where I’m making a convex hull for just a few points at once, but 10 - 50 million times. So need everything to be fast.
Can you elaborate?
Done! Still its not in Meshes.jl style but (if you are OK with this ) it might be very useful to somebody who needs to create a lot of small convex hulls.
For those not following the repository, this is the PR that @ig-or opened:
Thank you for the attempt @ig-or , I left some comments over there.
Update for those coming late here:
Please replace GeoTables.load by GeoIO.load. The function has been migrated to the GeoIO.jl package a long time ago.
I’m getting an error with the in operator. I think partition(::Any, ::Integer) in iterators.jl expects an Integer but is getting a Float64.
import GeoStats
import GeoIO
# file = "path/to/file" # download from "https://github.com/datasets/geo-countries/blob/master/data/countries.geojson"
geotable::GeoIO.GeoTable = GeoIO.load(file)
multipolygon = geotable.geometry[1]
# multipolygon isa GeoStats.Meshes.MultiPolygon
point = GeoIO.Meshes.Point(0, 0)
# point isa GeoStats.Meshes.Point
GeoStats.in(point, multipolygon) # throws error
GeoIO.in(point, multipolygon) # throws error
point in GeoStats.domain(geotable)[1] # throws error
The error
ERROR: MethodError: no method matching partition(::Meshes.Point2, ::Float64)
Closest candidates are:
partition(::Any, ::Integer)
@ Base iterators.jl:1269
Stacktrace:
[1] split_into_chunks(coll::Meshes.Point2, sz::Float64)
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:353
[2] _tcopy(xf::Transducers.Map{…}, ::Type{…}, reducible::Meshes.Point2, ::Transducers.SizeStable, ::Base.HasLength; basesize::Float64, kwargs::@Kwargs{})
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:347
[3] _tcopy(xf::Transducers.Map{…}, ::Type{…}, reducible::Meshes.Point2, ::Transducers.SizeStable, ::Base.HasLength)
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:345
[4] tcopy(xf::Transducers.Map{Meshes.var"#271#272"{Meshes.Box{…}}}, T::Type, reducible::Meshes.Point2; kwargs::@Kwargs{})
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:343
[5] tcopy(xf::Transducers.Map{Meshes.var"#271#272"{Meshes.Box{2, Float32}}}, T::Type, reducible::Meshes.Point2)
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:343
[6] tcollect(xf::Transducers.Map{Meshes.var"#271#272"{Meshes.Box{2, Float32}}}, reducible::Meshes.Point2; kwargs::@Kwargs{})
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:423
[7] tcollect(itr::Base.Generator{Meshes.Point2, Meshes.var"#271#272"{Meshes.Box{2, Float32}}}; kwargs::@Kwargs{})
@ Transducers ~/.julia/packages/Transducers/4xWio/src/reduce.jl:424
[8] sideof(points::Meshes.Point2, object::Meshes.Ring{2, Float32, CircularArrays.CircularVector{Meshes.Point2f, Vector{…}}})
@ Meshes ~/.julia/packages/Meshes/2WI3K/src/sideof.jl:79
[9] in(p::Meshes.Point2, poly::Meshes.PolyArea{2, Float32, Meshes.Ring{2, Float32, CircularArrays.CircularVector{…}}})
@ Meshes ~/.julia/packages/Meshes/2WI3K/src/predicates/in.jl:159
[10] (::Meshes.var"#251#252"{Meshes.Point2})(g::Meshes.PolyArea{2, Float32, Meshes.Ring{…}})
@ Meshes ~/.julia/packages/Meshes/2WI3K/src/predicates/in.jl:168
[11] _any
@ ./reduce.jl:1220 [inlined]
[12] any
@ ./reducedim.jl:1020 [inlined]
[13] in(p::Meshes.Point2, m::Meshes.MultiPolygon{2, Float32, Meshes.PolyArea{2, Float32, Meshes.Ring{…}}})
@ Meshes ~/.julia/packages/Meshes/2WI3K/src/predicates/in.jl:168
[14] top-level scope
Please open an issue with a MWE. It is hard to reproduce the error otherwise.
I’d be happy to open an issue. But I wanted to make sure first that I’m using the package as intended. At first it looks like the error is due to the distinction between Point vs Point2f
import GeoIO
import Meshes
# file = "path/to/file" # download from "https://github.com/datasets/geo-countries/blob/master/data/countries.geojson"
geometry = GeoIO.load(file).geometry # GeometrySet{2,Float32}
point32 = Meshes.Point2f(65, 35)
point64 = Meshes.Point(65, 35)
Meshes.in(point32, geometry[1]) # correctly false
Meshes.in(point32, geometry[2]) # correctly true
Meshes.in(point64, geometry[1]) # error
Meshes.in(point64, geometry[2]) # error
But that does not apply if I construct simple polygons myself:
p32 = Meshes.Point2f(0, 0)
p64 = Meshes.Point(0, 0)
polygon32 = Meshes.Ngon(Meshes.Point2f(-1, -1), Meshes.Point2f(1, -1), Meshes.Point2f(1, 1), Meshes.Point2f(-1, 1))
polygon64 = Meshes.Ngon(Meshes.Point(-1, -1), Meshes.Point(1, -1), Meshes.Point(1, 1), Meshes.Point(-1, 1))
p32 in polygon32 # correctly true
p32 in polygon64 # correctly true
p64 in polygon32 # correctly true
p64 in polygon64 # correctly true
If you think this is a bug, I will open an issue. I’m not sure I can reduce the MWE more than this though as I’m not familiar with the packages at all.
Our assumption in most algorithms is that all objects represent points with the same precision. If you mix float32 with float64 you will likely encounter errors. We can relax the dispatch in the future and promote to 64 bits, but that will need to wait the major CRS refactor that is ongoing.
Ok that makes sense. Thanks for your responses and for your work on the ecosystem.
1 Like
Does the inpolygon function still exist? I tried to use it today, but couldn’t find any trace of this function.
Try GeometryOps.jl, it has all the standard spatial predicates like within, contains etc that should do what you need:
3 Likes
Update regarding the Meshes.jl comments above. It now has a very efficient multi-threaded algorithm for point in polygon by Hao et al:
2 Likes