JuliaGeometry / Meshes.jl

Computational geometry in Julia
https://juliageometry.github.io/MeshesDocs/stable
Other
394 stars 85 forks source link

point in polygon errors when Point has LatLon as CRS #900

Closed disberd closed 3 months ago

disberd commented 3 months ago

Checking inclusion of a Point inside a Polygon with point in polygon throws an error if the point and the polygon are specified with LatLon as underlying CRS.

Here is a minimal reproducer:

using Meshes
# Generate PolyArea from set of coords
function create_polyarea(f, raw_coords)
    pts = map(raw_coords) do (;x,y)
        f(x,y)
    end
    r = Ring(pts)
    PolyArea(r)
end

# Function to generate normal Point
f_point(x,y) = Point(x,y)
# Function to generate LatLon point
f_latlon(lon, lat) = LatLon(lat, lon) |> Point

# Simple square coordinates
square_coords = map(45:90:360) do φ
    y,x = sincosd(φ)
    (;x,y)
end

# Tests
f_point(0,0) in create_polyarea(f_point, square_coords) # This returns true
f_latlon(0,0) in create_polyarea(f_latlon, square_coords) # This throws an error which is shown below

Picture of the error stack from Pluto: image

juliohm commented 3 months ago

As discussed on Zulip, we need to dispatch on Point{<:LatLon} and pick a valid algorithm from #615 for spherical/elliptical coordinates. In our first round of integration with CoordRefSystems.jl, we made sure that all algorithms for Cartesian coordinates worked as before. Now, in a second phase, we need to fix specific functions with other CRS.