Intersection of Ball and GeometrySet in LatLon crs

Specific Domains Geo
,
1

Hi, I’d like to find an intersection of a Ball with GeometrySet in LatLon coordinate reference system. However, when I construct my Ball b like

b = GeoStats.Ball(GeoStats.Point(LatLon(48.34, 18.45)), 10)

and ropes r as

38 GeometrySet ├─ Rope((lat: 49.7838°, lon: 18.2359°), ..., (lat: 49.7869°, lon: 18.2258°)) ├─ Rope((lat: 49.7869°, lon: 18.2258°), ..., (lat: 49.7946°, lon: 18.2327°)) ├─ Rope((lat: 49.7946°, lon: 18.2327°), ..., (lat: 49.795°, lon: 18.2337°)) ├─ Rope((lat: 49.795°, lon: 18.2337°), ..., (lat: 49.7973°, lon: 18.2488°)) ├─ Rope((lat: 49.7973°, lon: 18.2488°), ..., (lat: 49.7974°, lon: 18.2495°)) ⋮ ├─ Rope((lat: 49.8325°, lon: 18.1762°), ..., (lat: 49.8329°, lon: 18.176°)) ├─ Rope((lat: 49.8329°, lon: 18.176°), ..., (lat: 49.8351°, lon: 18.1754°)) ├─ Rope((lat: 49.8351°, lon: 18.1754°), ..., (lat: 49.8353°, lon: 18.1747°)) ├─ Rope((lat: 49.8353°, lon: 18.1747°), ..., (lat: 49.8345°, lon: 18.1697°)) └─ Rope((lat: 49.8345°, lon: 18.1697°), ..., (lat: 49.8341°, lon: 18.17°)),

then b ∩ r results in following error:
ERROR: StackOverflowError: Stacktrace: [1] intersection(f::Function, g₁::Ball{🌐, GeodeticLatLon{…}, Quantity{…}}, g₂::GeometrySet{🌐, GeodeticLatLon{…}, Geometry{…}}) @ Meshes C:\Users\oslejsek\.julia\packages\Meshes\K84Ms\src\intersections.jl:100--- the last 1 lines are repeated 79982 more times --- Some type information was truncated. Use show(err) to see complete types.

Different types of projection didn’t help and resulted in the same error. Then I tried to make 2 balls close to each other and apply the intersection function, which results in same error.

Is there any workaround around this?
PS: I am using GeoStats.jl v0.72.0 for this (but I guess the problem will be within Meshes.jl package).

2

Hi @stepanoslejsek , thank you for sharing the issue. Geodesic geometries usually require custom geometric algorithms, which we plan to implement in the first semester of 2025.

It would be great if you could reduce the example to a pair with a Ball and a single Rope. Your example is trying to compute intersection between a Ball and a GeometrySet of Ropes.

I’d first try to check if the intersection works with Cartesian coordinates. Then, try again with LatLon coordinates. If the error message shows up again, please feel free to open an issue on GitHub.

3

Sure, here is MWE for that:

using GeoStats

lats = [49.7868267, 49.7868428, 49.7868608, 49.7868831, 49.7869015, 49.7869172, 49.7869369]
lons = [18.2258339, 18.2258248, 18.2258177, 18.2258101, 18.225806, 18.225805, 18.2258069]

r = Rope([GeoStats.Point(LatLon(lat, lon)) for (lat,lon) in zip(lats,lons)])
# make center of the ball at one point of the rope
b = GeoStats.Ball(GeoStats.Point(LatLon(lats[5], lons[5])), 3)

intersection = b ∩ r

I tried Cartesian coordinates as well but the error is still persistent. I’ll open a github issue.

1 Like
4

Now I see the issue. You are trying to intersect a 2D ball wity a 1D rope. We don’t have a method for this pair yet, even for Cartesian coordinates and flat geometries.

1 Like
5

As a side note, notice that you don’t need to prepend all geometries with the GeoStats prefix unless you are loading two modules that export the same names. We usually import Makie backends with an explicit prefix instead, assuming that most of the geometric work is consistently handled by our stack:

import GLMakie as Mke

That way you can switch between Makie backends in a single place:

import CairoMakie as Mke

Regarding the issue above, we could start by adding an intersection method between Segment and Ball. The case with Rope could then be added with iteration over its segments.

The case with Segment and Ball can be as simple as checking a distance between the center of the Ball and the two vertices of the Segment. If any vertex is inside the Ball, you can then return a sub-Segment.

Please let us know if you plan to work on it. We can help with the PR review.

6

Not sure if I’ll manage to work on that, but I’ll try to manage it. Should I still open up an issue?

7

Please open an issue in the Meshes.jl repository. We will look into it.

1 Like