Operations on Polygons

Polygons can be manipulated with various different operations: (1) adding, which returns the union of polygons, (2) subtracting, which returns the difference of two polygons, (3) intersecting, which returns the regions both polygons contain, and (4) inverting, which takes a single polygon and returns a new polygon that encloses all the area within the original polygons bounding box that is not contained in the original.

using SpeciesDistributionToolkit
using CairoMakie

Adding Polygons

Let's start with the simplest version: adding polygons. Let's load the polygons for two US states, Texas and Oklahoma

texas = getpolygon(PolygonData(OpenStreetMap, Places); place = "Texas")
oklahoma = getpolygon(PolygonData(OpenStreetMap, Places); place = "Oklahoma")

FeatureCollection with 1 features, each with 0 properties

We can plot them to show that they are two adjacent states (that don't overlap):

Code for the figure

lines(texas)
lines!(oklahoma)

By default, these are both downloaded as FeatureCollections with a single feature. We can add them as FeatureCollections to yield a single polygon that consists of all regions enclosed by both of them.

tex_and_ok = add(texas, oklahoma)

Polygon(Geometry: POLYGON ((-103.00244140625 36.5450782775879,-103.0 ... 879)))

which we can then visualize:

Code for the figure

lines(tex_and_ok)

We can also do this with the single features directly

add(texas[1], oklahoma[1])

Polygon(Geometry: POLYGON ((-106.635765075684 31.8694763183594,-106. ... 594)))

or the polygons those features contain

add(texas[1].geometry, oklahoma[1].geometry)

Polygon(Geometry: POLYGON ((-106.635765075684 31.8694763183594,-106. ... 594)))

which yields the same result.

Alternatively, we can also just use the + operator, which does that same thing as add:

texas + oklahoma

Polygon(Geometry: POLYGON ((-103.00244140625 36.5450782775879,-103.0 ... 879)))

Sometimes we may have FeatureCollections with more than one feature. For example, lets combine Texas and Oklahoma FeatureCollections with vcat

vcat(texas, oklahoma)

FeatureCollection with 2 features, each with 0 properties

We'll also load another state to combine

louisiana = getpolygon(PolygonData(OpenStreetMap, Places); place = "Louisiana")

FeatureCollection with 1 features, each with 0 properties

By default, add unions all features within collections too, e.g. the following creates a single polygon with the entire area of all three states

tx_ok_la = texas + oklahoma + louisiana

Polygon(Geometry: POLYGON ((-91.6628646850586 29.3747406005859,-91.6 ... 859)))

which we can visualize:

Code for the figure

lines(tx_ok_la)

If we want to keep their separate boundaries, we should use vcat:

tx_ok_la_separate = vcat(texas, oklahoma, louisiana)

FeatureCollection with 3 features, each with 0 properties

And plot to show they keep their individual boundaries

fig-tx-ok-la-separate

lines(tx_ok_la_separate)

This covers the addition of polygons.

Subtracting polygons

We subtract polygons in a similar way. The method subtract(x, y) removes the area contained in y from x.

For example, let's remove the two largest metropolitan regions from Texas, first by loading them

dallas = getpolygon(PolygonData(OpenStreetMap, Places); place = "Dallas")
houston = getpolygon(PolygonData(OpenStreetMap, Places); place = "Houston")

FeatureCollection with 1 features, each with 0 properties

We can use subtract them iteratively

rural_texas = texas
rural_texas = subtract(rural_texas, dallas)
rural_texas = subtract(rural_texas, houston)

Polygon(Geometry: POLYGON ((-106.635757446289 31.8721027374268,-106. ... 664)))

and visualize:

Code for the figure

lines(rural_texas)

Or in a single line using the - operator

rural_texas = texas - dallas - houston

Polygon(Geometry: POLYGON ((-106.635757446289 31.8721027374268,-106. ... 664)))

which results in the same thing:

fig-rural-texas-sub-symbol

lines(rural_texas)

These operations can be chained to perform more complex clipping.

Intersecting a Polygon

We can also use the intersect method to extract the region where two polygons intersect.

As an example, let's consider trying to make a map of Texas, which different polygons corresponding to different ecoregions. We'll start by downloading the EPA North America 2nd level ecoregions

ecoregions = getpolygon(PolygonData(EPA, Ecoregions); level = 3)

FeatureCollection with 182 features, each with 3 properties

and then intersect them

tx_ecoregions = intersect(texas, ecoregions)

FeatureCollection with 12 features, each with 3 properties

and visualize the different ecoregions in texas

Code for the figure

poly(tx_ecoregions)

Finally, we can clip a series of polygons with a boundingbox:

bbox = (left=-100., right=-95.0, bottom=27.5, top=32.5)

tx_clip = clip(tx_ecoregions, bbox)

FeatureCollection with 8 features, each with 3 properties

fig-tx-clipped-by-bbox

lines(tx_clip)