Inside is a claim about the edges

Is this point inside that polygon is the second-most-asked question in spatial data, after which is nearest. It has a textbook answer — count the crossings of a ray, and an odd count means inside — and the answer is exact arithmetic on the stored coordinates.

The arithmetic is exact and the question is not well posed. A polygon is a list of vertices, and what joins them is not in the file.

What the file contains and what it does not #

A polygon in a spatial format is a sequence of coordinate pairs and a statement that the sequence is closed. Everything else is convention.

Two conventions are in use and they are not compatible:

• the edges are straight in the stored plane — join the vertices with straight lines in whatever coordinate system the file declares, which for a lon/lat file is the plate carrée;
• the edges are geodesics on the body — join the vertices along the ground, which is what somebody digitising a boundary from a map believed they were doing.

A straight segment is a claim about a plane measures the gap between those two readings for a single edge: the southern edge of the triangle above is 6,232 kilometres long and its two readings are 311 kilometres apart at the middle. This essay measures what that does to the interior.

The band, and its size #

Sampling the bounding box at 160 by 160 and asking both readings about every point gives the disputed set directly. For the triangle above: 3,590 of 25,600 test points disagree, and weighting each by its own area element gives 4.24 million square kilometres out of 14.40 — 29.5 per cent of the polygon.

Nearly a third. Not a boundary effect in any useful sense of the phrase: a third of the shape is in dispute between two readings of one file.

It is worth pausing on what that number is not. It is not an error bar on a measurement, because nothing here was measured — the vertices are exact and both readings compute their answers exactly. It is not a disagreement between two datasets, because there is one dataset. And it is not a defect in either algorithm, because both are correct implementations of what they are implementations of. It is a gap in what the data says, and the two answers are what fills the gap when a program has to produce one.

The size falls fast as the polygon shrinks, because it is governed by the edges’ own departure and that is quadratic in their length:

At the smallest of those the disputed band is still 2,200 square kilometres, which is a county. The share falls and the absolute area falls much faster, and neither reaches zero at any size a polygon is actually stored at.

Why a band rather than a line #

Two curves between the same two points bound a lens-shaped region, and the point-in-polygon answers differ exactly inside that lens. So the disagreement is two-dimensional — an area rather than a boundary uncertainty — and its width at the middle of an edge is the edge’s own departure.

That is why the share of the polygon is so large for a big triangle. The departure is 311 kilometres at the middle of a 6,232-kilometre edge, which is 5 per cent of the edge’s length; a lens 5 per cent as wide as it is long, on each of three edges, is a substantial fraction of a triangle.

For a polygon with many short edges the lenses are tiny and overlap nothing. The disputed area is set by the longest edge, not by the polygon’s size, which is the practical rule and the one a data producer can act on.

Densifying removes it, and is the only thing that does #

Inserting vertices along the intended curve makes the two readings agree, because a densified boundary is a polygon whose edges are short enough for the lens to vanish.

Measured: densifying the same triangle to about forty vertices per edge drops the disputed share from 29.2 per cent to 0.014 per cent — a factor of two thousand — and the residual is the sampling grid’s own resolution rather than any remaining geometry.

That is the control that makes the whole measurement a measurement. If the disagreement had survived densification it would be a bug in the point-in-polygon test rather than a fact about the polygon.

It is also the fix, and it is the only fix that works in a format that cannot express intent. A file cannot say “these edges are geodesics” in most schemas; it can say “here are more vertices”, and enough of them make the question moot.

A straight segment is a claim about a plane computes the vertex count a stated tolerance buys, and the law is a square root: a hundredfold tighter tolerance costs ten times the vertices. The same arithmetic applies here with the tolerance being the width of the band a producer is willing to leave in dispute.

The plane matters as well as the joining rule #

Everything above uses the plate carrée as the stored plane, because that is what a lon/lat file is. A polygon stored in a projected coordinate system has its own straight lines, and they are different straight lines.

So there are not two readings but as many as there are planes, and the pairwise disagreements between them are all of the same order. That is the same structure as a centroid belongs to a plane one operation along: a plane-dependent answer to a question about the ground, with no plane privileged by anything in the data.

Where in the polygon the disagreement is #

The band hugs the edges and its width varies along each one, so the risk of getting an answer wrong is not uniform over the shape.

At a vertex the two readings coincide exactly, because both pass through the stored point. At the middle of an edge the gap is widest. So the disputed region is three lenses, each pinched to nothing at both ends, and a point near a vertex is safe while one halfway along a long edge is not.

That has a consequence for testing. A test suite that checks containment at a polygon’s corners — which is what a person writing test cases naturally reaches for — will find the two readings agreeing perfectly and conclude there is no issue. The disagreement is at the midpoints of the longest edges, which is exactly where nobody puts a test point.

Which operations inherit it #

Every operation whose answer depends on the interior, which is most of them.

Containment queries — the case above.

Intersection and union. Two polygons’ overlap is a polygon whose boundary depends on the same convention, and its area therefore inherits the whole band. Two datasets from different producers, with different conventions, produce sliver polygons along every long shared edge.

Area. Computing an area needs a surface measures what the plane does to an area computed by shoelace, and this adds a second term: even in the right plane, the area depends on which curve the edges are.

Rasterisation. A grid built by asking “is this cell’s centre inside” is exactly this question repeated, so the raster’s edge is the band’s edge, and two rasterisations of the same polygon differ in whole rows near a long edge.

What a projection does not fix #

It is worth being explicit that no choice of projection removes this, because the instinct on a site about projections is to look for one that does.

Reprojecting the polygon reprojects its vertices. The edges are still whatever the reader takes them to be, and the ambiguity travels with the data unchanged — a polygon reprojected from lon/lat into an equal-area plane has the same three vertices and the same undeclared edges, so the same two readings exist there and disagree by the same ground.

What a projection does change is which reading is cheap. On the gnomonic projection a geodesic edge is a straight line, so a polygon whose edges are geodesics can be tested by ordinary plane arithmetic there — which makes the gnomonic the natural plane for exact spherical containment over a hemisphere, and is one more use for the gnomonic companion’s single exact property.

That is a real answer to the problem for one convention and it does not scale: the gnomonic shows less than half the sphere, so a polygon spanning more than that cannot be tested in one plane at all, and the convention has to be honoured by densification instead.

The two readings of a real boundary #

Long straight edges are not a contrivance. They are how a great many real boundaries are defined: a treaty line along a parallel, a survey line along a meridian, a maritime limit drawn as a straight segment between two coordinates, a management zone drawn on a chart with a ruler.

For those the question this essay asks is not academic and the answer is legal rather than geometric. A boundary defined as “the parallel of 49° north” is a curve on the ground — a small circle, not a geodesic and not a straight line in any projected plane — and none of the three readings above is it. A boundary defined as “the straight line between these two points” needs a plane named, and the treaty rarely names one.

The geometry’s contribution to that is small and specific: it can say how far apart the readings are, so that a reader knows whether the ambiguity matters. A segment along the 50th parallel 143 kilometres long has its two readings 478 metres apart at the middle; one 1,425 kilometres long has them 48 kilometres apart. The quadratic law again, and the crossover between “a surveying detail” and “a territorial dispute” sits somewhere in between.

The convention actually used #

There is no universal answer and the situation is worse than merely unstandardised, because different parts of one stack differ.

Most desktop and library software treats edges as straight in the stored plane, because that is what a plane sweep is. Some spatial databases offer a geography type whose edges are geodesics and a geometry type whose edges are not, in the same database, with implicit conversion between them. Some formats’ specifications state a convention and most implementations of those formats do not check it.

The site takes no view on which is right — both are defensible and the choice is a matter for whoever defines the data — and states the two things geometry can settle: the two readings differ by an area, the area is set by the longest edge, and densification makes the question go away.

What a producer should do #

Three things, in order of how much they cost.

Densify long edges before publishing. A tolerance in metres and the square-root law give the vertex count, and the cost is a few times the file size for a boundary whose longest edge is a continent.

State the convention in the metadata even where the schema has no field for it. A reader who knows which reading was meant can reproduce it; one who does not cannot.

Do not mix. A dataset assembled from producers with different conventions has slivers along every long shared edge, and the slivers are not detectable after the fact — they look exactly like a genuine disagreement about where a boundary is. That is the same class of silent inheritance reprojecting a raster invents values describes for imagery: an operation that had to make something up, with nothing in the output saying it did.

What it costs to do properly #

The exact operation for the geodesic reading is not expensive, which is worth saying because “do it properly” usually means “do it slowly”.

Testing a point against a spherical polygon is the same crossing count with spherical arithmetic: for each edge, decide which side of the great circle through its endpoints the point lies on, which is one dot product with the edge’s pole. That is three multiplications and two additions per edge, against two multiplications and a division for the planar test.

So the exact version costs roughly what the approximate one costs, and the reason nobody uses it is not performance. It is that the planar test is what every library shipped, and it returns an answer that is right nearly everywhere.

That is the same conclusion nearest is a question about the metric reaches for the neighbouring operation, with one difference worth noting: there the exact computation is genuinely more expensive — an ellipsoidal geodesic is an iteration — and the answer is a filter-then-decide structure with a bounded error. Here there is no such excuse.

The band against the file’s own precision #

There is one comparison that decides whether a producer needs to care, and it can be made from the two numbers a file already carries: how long its edges are, and how many decimal places it stores.

The band’s half-width is the sagitta of the edge — the departure of the chord from the arc — which is $\kappa L^2/8$ for an edge of length $L$ on a sphere of radius $R$, with $\kappa = 1/R$. Setting that against a stored precision $p$ gives a crossover length:

$$L^\star = \sqrt{8 p R}$$

At seven decimal places of degrees — the precision of a great many published files, and about a centimetre on the ground — that is 714 metres.

The reading is sharp on both sides of it.

Below 714 metres, the ambiguity is smaller than the file’s own rounding. A dense boundary with vertices a few hundred metres apart is unambiguous at the precision it is stored to: the two interpretations of every edge differ by less than the last digit, so no consumer can tell them apart and none needs to be told.

Above it, the ambiguity grows as the square while the precision stays fixed, and the gap opens fast. A ten-kilometre edge has a two-metre band, which is two hundred times the stored precision; a hundred-kilometre edge has a 196-metre band, twenty thousand times it.

So a file with kilometre-scale edges and centimetre-scale coordinates is stating a position to a precision it cannot mean, by three to four orders of magnitude — and it is doing so in the digits, where every consumer reads it as an assurance. That is the defect in its most reportable form, and it needs no test point and no library: the ratio of $\kappa L^2/8$ to the stored precision is computable from the file alone.

It also gives the producer’s remedy a number rather than an adjective. Densify is advice; densify until no edge exceeds 714 metres is a specification, it follows from the precision the producer already chose to publish at, and it is the length at which the two readings of the file stop being two.

Where this ladder goes #

The dataset ladder has established what a coordinate refers to, what a degree is not, what an area needs, what a straight segment claims, where the antimeridian cuts, which neighbour is nearest, what a reprojection invents, and which plane a centroid belongs to.

This adds the interior: the operation that decides membership, and the discovery that the answer is a band rather than a point because a file records vertices and not edges.

What that leaves is the shape of the whole field’s argument, which the operation decides the coordinate system states: no stored system makes every operation exact, and the reason is the theorem this site opened with. This essay adds a corollary the theorem does not cover — some operations are not made exact by any system, because the ambiguity is in what the data means rather than in the plane it is stored in.