Tutorial

Point inclusion algorithm

This package was born from the need of ensuring a simple way of optimizing the performance of checking whether a point on the Earth's surface is located within a given region (e.g. a country, a satellite beam, or any other area that can be represented by a polygon in latitude/longitude coordinates).

We need to define multiple types of such regions in our downstream packages (e.g. CountriesBorders.jl) and we want to minimize code duplication as much as possible.

The traditional p in region is well defined in Meshes.jl but is suboptimal in terms of performance when one wants to check a lot of points (e.g. thousands) over a domain spanning multiple polygons (e.g. the domain of all polygons associated to all countries). The complexity of checking for point inclusion in a polygon increases with the number of points in the polygon, and when you have lots of points and lots of polygons, most of your point/polygon pair will return false.

The approach to speed up computations taken in this package is to create custom geometries that store not only the polyareas but also the bounding boxes (i.e. a Box from Meshes.jl) associated to each polyarea. The modified inclusion algorithm then simply prefilters point/polygon pairs by checking first inclusion in the bounding box (which is extremely fast) and falling back to checkin inclusion in the polygon only if p in bbox is true. This has shown to provide speed ups of 20x-40x in tests when implemented for the CountriesBorders.jl downstream package (see this PR)

FastInGeometry Interface

To simplify exploiting this fast algorithm for custom geometries defined in downstream packages, this package defines an interface that is adhered to when a type satisfies the three conditions below:

• It subtypes the FastInGeometry abstract type defined and exported by this package
• It has a valid method for the interface function polyareas returning an iterable of PolyAreas.
• It has a valid method for the interface function bboxes returning an iterable of Boxs

In reality, the 2nd and 3rd conditions can also be met by simply containing a field which subtypes GeoBorders (more details in the next section) or defining a custom method for the interface function geoborders that returns an instance of type GeoBorders.

As an example, see the following code snippet defining a simple geometry that satisfies the FastInGeometry interface by resorting to the simplest approach of having a field subtyping GeoBorders

using GeoBasics

# Define a custom type which subtypes 
struct SimpleGeometry <: FastInGeometry{Float64}
    name::String
    borders::GeoBorders{Float64}
end

GeoBorders

This package defines a single concrete subtype satisfying the FastInGeometry interface with the GeoBorders type. This is intended to represent the borders of arbitrary regions, and provides an easy way for custom regions to adhere to the FastInGeometry interface (as explained above).

Additionally, when provided with geometries that are plain Multi, PolyArea or Box from Meshes (and not other geometries satysfiying the FastInGeometry interface), the constructor for GeoBorders automatically tries fixing issues with polygons crossing the antimeridian line. It does so by implementing the algorithm of the antimeridian python library that also provides a simplified explanation of the algorithm in its documentation.

To show the problem of with polygons crossing the antimeridian, consider the following complex polygon that contains holes and crosses the antimeridian multiple times

using GeoBasics
using GeoBasics.Meshes
using GeoBasics.GeoPlottingHelpers
using PlotlyBase
using PlotlyDocumenter

complex_s_poly = let
    f = to_cartesian_point
    outer = map(f, [ # Exterior part of the polygon, crossing the antimeridian multiple times
        (160,30),
        (-160,30),
        (-160,0),
        (170, 0),
        (170, -10),
        (-160, -10),
        (-160, -20),
        (160, -20),
        (160, 10),
        (-170, 10),
        (-170, 20),
        (160, 20),
    ]) |> Ring
    inner_west = map(f, [ # Hole in the western hemisphere, in the lower right of the polygon
        (-162, -12),
        (-162, -18),
        (-168, -18),
        (-168, -12)
    ]) |> Ring
    inner_east = map(f, [ # Hole in the eastern hemisphere, in the upper left of the polygon
        (162, 22),
        (162, 28),
        (168, 28),
        (168, 22)
    ]) |> Ring
    inner_both = map(f, [ # Hole crossing the antimeridian, in the middle of the polygon
        (170, 2),
        (170, 8),
        (-170, 8),
        (-170, 2)
    ]) |> Ring
    PolyArea([outer, inner_west, inner_east, inner_both])
end

plt = with_settings(:OVERSAMPLE_LINES => :SHORT) do # This makes sure that the plot uses the shortest line between points
    data = geo_plotly_trace(complex_s_poly)
    Plot(data, Layout(;
      geo = attr(;
        center_lon = 180,
        lonaxis_range = [0, 360]
      )
    ))
end

In reality, the polygon looks fine because we forced the plot to use shortest line between points for plotting (thus crossing the antimeridian). The actual polygon defined with those coordinates represents instead the following degenerate region:

plt = with_settings(:OVERSAMPLE_LINES => :NORMAL) do # Actually draw lines from endpoints without using shorted line
    data = geo_plotly_trace(complex_s_poly)
    Plot(data)
end

which can also be verified by checking point inclusion of one point that shouldn't be inside but is, and of one that should be in but isn't:

# We use cartesian point as point inclusion in LatLon is not always defined in Meshes
should_not = to_cartesian_point(LatLon(25, 0)) # This is in africa and shouldn't be in the intended polygon
should_be = to_cartesian_point(LatLon(25, 180)) # This is in the ocean inside the polygon and outside it's holes

map(in(complex_s_poly), (;should_be, should_not))

(should_be = false, should_not = true)

By feeding this polygon into the constructor of GeoBorders, the antimeridian crossin is handled by splitting the input polygon into 4 subpolygons whenever a crossing of the antimeridian is encountered:

gb = GeoBorders(complex_s_poly)
plt = with_settings(:OVERSAMPLE_LINES => :NORMAL) do # Actually draw lines from endpoints without using shorted line, showing that GeoBorders handle this correctly
    data = geo_plotly_trace(gb)
    Plot(data, Layout(;
      geo = attr(;
        center_lon = 180,
        lonaxis_range = [0, 360]
      )
    ))
end

And checking again for point inclusion we find the expected behavior (including a point in the hole not being considered inside).

inside_hole = LatLon(25, 165) # Notice that this does not need to be a Point

map(in(gb), (; should_be, should_not, inside_hole))

(should_be = true, should_not = false, inside_hole = false)

As seen from the last example, there is one additional advantage to satisfying the FastInGeometry interface. Inclusion in a FastInRegion does not need to use points which are of Point type with the exact same CRS as the geometry (as in plain Meshes) but can be:

• a Point with either LatLon{WGS84Latest} or Cartesian2D{WGS84Latest} CRS,
• a plain LatLon{WGS84Latest} coordinate.

Override antimeridian fix

The procedure in the GeoBorders constructor that fixes the antimeridian crossing decides if there is one by simply checking if any segment of a polygon spans more than 180° degrees of longitude.

In some cases, this has the unintended consequence of modifying the intended polygon. Consider for example a rectangular polygon going from -100 to +100 longitude and from -20 to + 20 latitude. Creating a GeoBorders with this input polygon will uncorrectly split it:

large_rectangle = let
  f = to_latlon_point(Float64)
  PolyArea(map(f, [
    LatLon(-20, -100),
    LatLon(-20, 100),
    LatLon(20, 100),
    LatLon(20, -100),
  ]))
end

plt = with_settings(:OVERSAMPLE_LINES => :NORMAL) do # Actually draw lines from endpoints without using shorted line
    data = geo_plotly_trace(GeoBorders(large_rectangle))
    Plot(data)
end

In these cases, it is possible to override this behavior by using the fix_antimeridian_crossing keyword argument:

plt = with_settings(:OVERSAMPLE_LINES => :NORMAL) do # Actually draw lines from endpoints without using shorted line
    data = geo_plotly_trace(GeoBorders(large_rectangle; fix_antimeridian_crossing = false))
    Plot(data)
end

FastInDomain Interface

As alternative to subtyping FastInGeometry, custom types can also participate to the fast inclusion algorithms if they subtype FastInDomain.

In this case, the custom subtype must represent a domain containing geometries implementing the FastInGeometry interface and must implement have valid methods for the following two functions from Meshes.jl:

• Meshes.nelements(custom_domain): This should return the number of geometries within the domain
• Meshes.element(custom_domain, ind): This should return the ind-th geometry from the domain

See an example implementation of a custom domain with elements of a custom geometry satisfying the interface

struct SimpleGeometry <: FastInGeometry{Float64}
    name::String
    borders::GeoBorders{Float64}
end

function SimpleGeometry(geoms; name = "SimpleGeometry")
    borders = GeoBorders{Float64}(geoms)
    return SimpleGeometry(name, borders)
end

struct SimpleDomain <: FastInDomain{Float64}
    name::String
    geoms::Vector{SimpleGeometry}
end

Meshes.nelements(sd::SimpleDomain) = length(sd.geoms)
Meshes.element(sd::SimpleDomain, ind::Int) = sd.geoms[ind]