Create array of hexagons, place in index, and search by point

Specific Domains Geo
1

I’ve been looking at JuliaGeo to how I might do the following: 1) create a tiled array of hexagons, 2) put them in a spatial index, and 3) be able to find which hexagon intersects a point. I have written this code in Boost Geometry and would like to see how it looks in Julia.

I assume I can use a LineString or LinearRing to build the hexagons. Can you give an example of constructing these that uses pure Julia (ie not pointers into eg GDAL)?

How would I then go about bulk-loading these into an rstar tree and doing point queries?

Thanks.

2

@thk686 if you are interested in a pure Julia solution, I recommend checking the GeoStats.jl stack instead, particularly the Meshes.jl submodule.

You can create an hexagon given a list of vertices in counter-clockwise orientation:

julia> using Meshes

julia> points = [(cos(θ),sin(θ)) for θ in range(0, 2π-π/3, length=6)]
6-element Vector{Tuple{Float64, Float64}}:
 (1.0, 0.0)
 (0.4999999999999999, 0.8660254037844387)
 (-0.5000000000000002, 0.8660254037844386)
 (-1.0, 1.2246467991473532e-16)
 (-0.4999999999999996, -0.8660254037844388)
 (0.5000000000000001, -0.8660254037844386)

julia> hexagon = Hexagon(points)
Hexagon{2,Float64}
  └─Point(1.0, 0.0)
  └─Point(0.4999999999999999, 0.8660254037844387)
  └─Point(-0.5000000000000002, 0.8660254037844386)
  └─Point(-1.0, 1.2246467991473532e-16)
  └─Point(-0.4999999999999996, -0.8660254037844388)
  └─Point(0.5000000000000001, -0.8660254037844386)

and can visualize for debugging purposes:

using MeshViz
import GLMakie as Mke

viz(hexagon)

image
image782×586 10.3 KB

If you are interested in finding out which hexagon is closest to a given point, you can place the hexagons inside a Collection, which is a “domain” in our jargon:

gset = Collection([hexagon1, hexagon2, ...])

Then use the KNearestSearch method to query points. Keep in mind that we are actively working on adding new methods, geometries and domain types. Contributions are very welcome.

3

Sweet. Very nice. I will try this. I can use nn search to find the enclosing polygon when the point falls inside the area of interest, but more generally, I will also need an intersects predicate to know when the point falls outside any polygon. Is that planned or should I use eg GEOS for that part? I’d be happy to help out if I can.

1 Like
4

We already have intersection algorithms, and algorithms to check the relative position of a point in a geometry. Check

https://juliageometry.github.io/Meshes.jl/stable/algorithms/intersection.html

You can also use

point ∈ hexagon

to check if a point is in a polygon, and also

sideof(point, boundary(hexagon))

to check if it is :INSIDE, :OUTSIDE or :ON the boundary.

1 Like
5

Since this question concerns hexagons, floor and parity operations on the transformed point coordinates should identify the hexagon it belongs to without forming an array of hexagons and intersecting with each one as a point-to-polygon query.

The hexagon tiles consist of triangle tiles whose vertices form a grid with regular unit steps and a 60-degree angle instead of a right angle at the origin. That is how I would do it.

Welcome to the Julia language discourse!

2 Likes
6

Agree 100% @pitsianis , we could exploit the regularity of the tessellation to perform the queries. Unfortunately, we don’t have a HexagonGrid domain implemented yet. Or mechanisms to operate on these tilings. People have shared comments about this on Zulip, and I welcome everyone to join our forum over there.

7

Thanks for the comment. Yes, I know I could implement this as a type of grid. I want the code to be generalizable to any set of polygons, so it make sense to have a generalized spatial index.

8

My function generates a vector of hexagons. Is there a way to plot all at once? Also, it seems a Collection takes an array. Can I initialize a Collection from a vector?

9

Figured it out… [myvec[i] for i in eachindex(myvec)]

10

You can simply pass the vector of geometries to the viz function and it will recognize each individual geometry. The Collection takes a AbstractVector of geometries, you can Collection(collect(list)) if your list is not already a Vector.

11

What is the type of myvec? You can check if myvec isa AbstractVector. If it is not, you can collect(myvec), which is basically equivalent to the list comprehension syntax.

12

Thanks. Getting used to the type system. I realized after posting that it was a type issue and modified the function to return a concrete type.