Contour Models for Physical Boundaries Enclosing Star-Shaped and Approximately Star-Shaped Polygons

Abstract

Boundaries on spatial fields divide regions with particular features from surrounding background areas. Methods to identify boundary lines from interpolated spatial fields are well established. Less attention has been paid to how to model sequences of connected spatial points. Such models are needed for physical boundaries. For example, in the Arctic ocean, large contiguous areas are covered by sea ice, or frozen ocean water. We define the ice edge contour as the ordered sequences of spatial points that connect to form a line around set(s) of contiguous grid boxes with sea ice present. Polar scientists need to describe how this contiguous area behaves in present and historical data and under future climate change scenarios. We introduce the Gaussian Star-shaped Contour Model (GSCM) for modelling boundaries represented as connected sequences of spatial points such as the sea ice edge. GSCMs generate sequences of spatial points via generating sets of distances in various directions from a fixed starting point. The GSCM can be applied to contours that enclose regions that are star-shaped polygons or approximately star-shaped polygons. Metrics are introduced to assess the extent to which a polygon deviates from star-shapedness. Simulation studies illustrate the performance of the GSCM in different situations.

1 INTRODUCTION

Boundaries that enclose regions are often subjects of scientific interest. Contour lines divide a contiguous region with some defining feature(s) from surrounding background areas. This paper introduces the Gaussian Star-shaped Contour Model (GSCM) for the distribution of such contours. The GSCM is designed for modelling contours that enclose regions that are star-shaped polygons (Definition 4) or approximately star-shaped polygons.

The GSCM is motivated by the need for models appropriate for the sea ice edge. In the Arctic ocean, large areas are covered by sea ice, or frozen ocean water. The sea ice edge contour forms a boundary between the area covered by sea ice and the surrounding open water. Data on the location of sea ice is provided by remotely sensed gridded products that indicate if a grid box is ice-covered. Since most grid boxes with sea ice are in one contiguous region, this data can be converted to binary values that indicate if each grid box is inside or outside the main ice-covered area. The ice edge contour is then the set of points that connect to form the boundary between the grid boxes inside and outside the region. Thus the ice edge contour is modelled as a collection of ordered, connected spatial points in the two-dimensional (2D) plane.

Figure 1 shows a sample sea ice edge contour. Polar scientists are interested in the variability of ice edge contours and the extent to which ice edge contours change over time. Predictions of the ice edge contour are also needed weeks to months in advance for maritime planning. Distributions of ice edge contours are inferred from observing multiple ice edge contours, such as those observed at different times.

Previous research has developed contour and boundary models for other data types. Analysis has often focused on inferring a single boundary from observations on a spatial field. Research on exceedance levels has developed methods to infer contours that describe where a property goes above some level (Bolin & Lindgren, 2015; French & Hoeting, 2016; French & Sain, 2013). Wombling methods find contours by identifying curvilinear gradients (Banerjee & Gelfand, 2006; Womble, 1951). Statistical shape analysis (Dryden & Mardia, 2016; Srivastava & Klassen, 2016) provides tools to model boundaries corresponding to particular objects with discernible features. None of these methods infer distributions of contours from multiple observed contours like the sea ice edge data. Similarly, remote sensing research on sea ice has focused on how to infer the location of the sea ice edge at individual time points given information from active and passive remote sensing. So, the focus is again on how to infer a specific boundary from data. The GSCM addresses a different question: how to model the distribution of sea ice edges given samples of fully observed contours represented as connected sequences of points.

We first develop and assess the GSCM for modelling boundaries enclosing star-shaped and approximately star-shaped polygons. We then apply GSCMs to solve our motivating problem of how to model the sea ice edge. Section 2 defines contours and how to represent them. Section 3 introduces the GSCM for modelling contours enclosing star-shaped polygons and discusses model fitting. Section 4 introduces a metric for assessing coverage of star-shaped contours and Section 5 presents results from simulation studies. Section 6 extends the GSCM to contours enclosing approximately star-shaped polygons. Section 7 applies the GSCM to sea ice edge contours. Section 8 concludes the paper with discussion, including Section 8.1 that examines how the GSCM compares to other contour and boundary methods.

2 CONTOUR DEFINITIONS

Our focus is on modelling contours that act as the boundary between a region that has some feature(s) and the surrounding background region, that is, ice-covered areas versus open water. There are multiple ways such contours could be defined. In this section, we give two representations for these contours that will be used as a basis for subsequent modelling and assessment.

2.1 Point-sequence representation

Contours and the regions they enclose can be described using connected sequences of points. We refer to this description of a contour as a point-sequence representation. We define the following concepts.

Definition 1

Contour point sequence, $S$⁠: An ordered set of spatial points $S=(s1,…,sn)$⁠, with $n>2$⁠, where each $si$ consists of the x-y coordinates of a spatial location.

Definition 2

Contour line, $S‾$⁠: The connected line formed by connecting $si$ to $si+1$ for $i=1,…,n−1$ and connecting $sn$ to $s1$⁠.

Definition 3

Enclosed polygon, $S_$⁠: The polygon formed by the interior of the contour line $S‾$⁠.

The left panel of Figure 2 illustrates these definitions for a contour described by a point-sequence representation. The main advantage of the point-sequence representation is its flexibility. Any contour enclosing a polygon can be represented exactly with a sequence of spatial points, $S$⁠. Also, the level of detail represented can be increased simply by increasing the number of points.

A binary grid indicating membership in the region of interest can be converted to the point-sequence representation and modelled accordingly. The contour $S$ is made up of corner points of grid boxes that touch the outside of the region on one side and the inside of the region on the other. The points of $S$ are ordered to align with the order in which they would be touched if one were to trace around the boundary. Where and in what direction to start tracing around the boundary is arbitrary. These choices only determine the indexing of the points in $S$⁠, not the line, $S‾$⁠, or enclosed polygon, $S_$⁠. The point-sequence representation can be used with any grid resolution, though finer grids will require more points in $S$⁠. For sea ice, the binary grid is what is observed.

2.1.1 Notation

We need to distinguish between points, lines and polygons. An ordered sequence of spatial points will be denoted by a boldface letter, such as $S$⁠. A line formed by connecting points will be denoted with an overline, such as $S‾$⁠. A line segment will be denoted by an overline over two letters that represent the start and end points of the segment, such as $CD‾$⁠. The polygon enclosed by a line will be denoted by an underline, such as $S_$⁠.

2.1.2 Contours with fractal characteristics

We acknowledge that the point-sequence representation does not directly account for contours whose true nature is fractal. As such, representing fractal contours as connected sequences of points may be an approximation. In contours of sea ice and other physical world examples, as the spatial scale of observations increases, the level of detail of the contour also increases (Mandelbrot, 1967). With each increase in spatial resolution additional line segments are needed to describe the increased detail. In other words, the length of the contour increases each time the spatial resolution increases. This fractal nature of some contours means that these contours can never be fully expressed with a finite, ordered set of spatial points. These contours' true lengths are infinite.

For the applications of interest, however, limits exist on the precision of measurements and the relevant scale of scientific interest. While the fractal or Hausdorff dimension can be estimated (Gneiting et al., 2012), the level of detail of the boundary that will be measured or needed will rarely be a fractal. So, making the simplifying assumption that the contour can be defined by a sequence of spatial points is reasonable for modelling the sea ice edge contour. Additional discussion of contours as fractals is given in Section 8.2.

2.2 Star-shaped representation

Point-sequence representations are natural and describe contours accurately. However, point-sequence representations are ill-suited for describing multiple contours and their distributional behaviour. Contours differ in length, so two points with the same index on two different point-sequence representations are not likely to be in the same physical location. Comparing spatially dependent features and inferring distributions is therefore difficult. We build an alternate star-shaped representation that avoids the weaknesses of point-sequence representations. The star-shaped representation is appropriate for contours that enclose star-shaped polygons or approximately star-shaped polygons. The general idea of star-shaped representations is to describe contours based on the length of the lines extending to the contour from a fixed point in different directions rather than as a sequence of spatial points as in a point-sequence representation. Figure 2 illustrates how a point-sequence representation (left) and a star-shaped representation (right) differ.

Before defining the star-shaped representation, we review the standard definitions of a star-shaped polygon and its kernel (Preparata & Shamos, 1985, p. 18).

Definition 4

Star-shaped polygon: A polygon $P_$ is star-shaped if there exists a point $D$ within $P_$ such that the line segment $Dp‾$ is fully contained within $P_$ for all points $p$ on line $P‾$⁠.

All convex polygons are star-shaped, but the set of star-shaped polygons is substantially larger. Figure 3 shows nine example star-shaped polygons. These examples highlight the variety of polygons that can be star-shaped.

Definition 5

Kernel of a star-shaped polygon, $𝒦(P_)$⁠: The set of point(s) that satisfy the criterion for $D$ in Definition 4 is referred to as the kernel of the polygon, $𝒦(P_)$⁠.

Convex polygons are the subset of star-shaped polygons such that $𝒦(P_)=P_$⁠.

For any star-shaped polygon, $S_$⁠, lines can be drawn from some point, $C$⁠, in the kernel of $S_$ to all points on a contour, $S‾$⁠. Assume for the moment that the contour point sequence $S$ is unknown, but that the location of $C$ is known along with the lengths and directions of the lines from $C$ to $S$⁠. Then $S$ could be derived from this information with trigonometry. Taking inspiration from this fact, we develop the star-shaped representation. First we define a line set:

Definition 6

Line set, $ℒ(C,θ)$⁠: A set of $p>2$ lines, $ℒ=(ℓ1,…,ℓp)$⁠, extending infinitely outward from some starting point $C=(Cx,Cy)$ at $p$ unique angles, $θ=(θ1,…,θp)$⁠, where $Cx$ and $Cy$ are x- and y-coordinates, respectively.

We also define a set of spatial points that produce a star-shaped polygon when connected in order:

Definition 7

Star-shaped point set, $V(C,θ,y)$⁠: A set of $p>2$ spatial points $V=(v1,…,vp)$ such that

$vi=(Cx+yicos(θi),Cy+yisin(θi)),$

(1)

where $y=(y1,…,yp)$ is a set of $p$ distances, $C=(Cx,Cy)$ is a spatial point, and $θ=(θ1,…,θp)$ is a set of $p$ unique angles.

A star-shaped point set can be used to represent a contour when the distances are selected systematically:

Definition 8

Star-shaped representation, $V˜(C,θ,S$⁠): Let $S_$ be a star-shaped polygon and $C∈𝒦(S_)$ be a starting point. Then, the star-shaped representation of the contour $S$⁠, denoted by $V˜(C,θ,S)$⁠, is the star-shaped point set, $V(C,θ,y)$⁠, where $y=(y1,…,yp)$ is the set of distances from $C$ to the intersection point of the contour line $S‾$ and each line $ℓi$ in the line set $ℒ(C,θ)$⁠.

The right panel of Figure 2 shows the components of the star-shaped representation for a sample contour. Let $V‾$ refer to the contour line formed by connecting $vi$ to $vi+1$ for $i=1,…,p−1$ and $vp$ to $v1$ for $p>2$⁠. Let $V_$ refer to the polygon contained within $V‾$⁠.

Theorem 1

Let$θ=(θ1,…,θp)$with$θi<θi+1$and$θi∈(0,2π)$for all$i$⁠. For a star-shaped polygon$S_$there exist$θ$and$y$such that$V˜(C,θ,y)=S$for any$C∈𝒦(S_)$⁠. (Proof in Appendix A.)

Corollary 1

Let$ℓθ$denote the line that extends infinitely outward from$C$at angle$θ∈(0,2π)$and that intersects$S‾$⁠. For any$θ$⁠, the line$ℓθ$is distinct, that is,$ℓθ≠ℓθ′$for any$θ,θ′$such that$θ≠θ′$⁠. (Proof in Appendix A.)

The star-shaped representation allows for finding how contours differ and what the variability of the contours is in different spatial areas. For example, consider two contours $Sk$ and $Sℓ$ described with star-shaped representations, $V˜(C,θ,Sk)$ and $V˜(C,θ,Sℓ)$⁠, for common line set $ℒ(C,θ)$⁠. To find how much further one contour extends in any direction, simply find the difference between $yi,k$ and $yi,ℓ$ where $yi,k$ and $yi,ℓ$ are the distances from $C$ to contours $Sk$ and $Sℓ$ along a line extending at angle $θi$⁠. The variability of the contours along any line $ℓi∈ℒ(C,θ)$ is estimated from the variability of the corresponding $yi$ values in the contours' star-shaped representations.

For a contour enclosing a star-shaped polygon, the star-shaped representation is identical to the point-sequence representation when $p=n$⁠, $C∈𝒦(P_)$⁠, and $θi$ aligns with the direction of the line segments $Csi‾$ for all $i$⁠. When these conditions are met, the points $V$ are the same for any choice of starting point $C$ within the kernel of $S_$⁠. However, the angles $θ$ and lengths $y$ will differ depending on $C$⁠.

3 STAR-SHAPED CONTOUR MODEL

3.1 General model

We now propose the star-shaped contour model for generating contours that enclose star-shaped polygons. In later sections, we will build upon this idealised model for star-shaped polygons to develop a model appropriate for ice edge contours that are only approximately star-shaped.

Definition 9

Star-shaped contour model, ${C,θ,Γ}$⁠: Let $C=(Cx,Cy)$ be a fixed starting point, let $θ=(θ1,…,θp)$ be a fixed set of $p>2$ unique angles, and let $Γ$ be a probability distribution from which a set of values $y=(y1,…,yp)>0$ can be drawn. This set of parameters form a star-shaped probability model if each drawn set $y$ can be used to form a corresponding star-shaped points set, $V(C,θ,y)$⁠, as given in Definition 7.

We now consider a distribution that is appropriate in many circumstances. We assume that $y$ follows a Gaussian distribution,

where $μ$ is a mean vector and the parameter $∑$ is a positive-definite covariance function. We further assume that $μ$ and $∑$ are such that mass on non-positive $y$ is negligible. (In practice, if in a small proportion of cases, a generated $yi$ is non-positive, its value can be set to some small $η>0$⁠.) We call this model the GSCM. The GSCM can be seen as a finite approximation to the planar version of Gaussian Random Particles proposed in Hansen et al. (2015). GSCMs can produce a fairly flexible set of contours. The first three rows in Figure 3 illustrate some types of contours GSCMs can produce.

Because of how $y$ is constructed, reasonable $μ$ and $∑$ that align with typical observable $y$ values will avoid substantial non-positive $y$⁠. The values $y$ represent distances from a starting point $C$⁠, so are automatically non-negative. Some points in the kernel of the polygon will typically be centrally located points. Lengths close to zero can be avoided by using one of these points for $C$⁠.

Covariance matrices, $∑$⁠, are based on the structure of the lines in the line set. The correlation structure for the set of distances, $y$⁠, is a function of the angles, $θ$⁠, of the lines in the line set. Covariances based on angles are complicated by the fact that 0 and $2π$ represent the same angle. So, the difference between two angles does not necessarily correspond to how far apart the angles actually are. Specialised covariance functions have been derived that remain valid when distances are indexed by angle (Gneiting, 2013). Denote the angle between $θi$ and $θj$ by $d(θi,θj)∈[0,π]$⁠.

Typically, the correlation between $yi$ and $yj$ will decrease as $d(θi,θj)$ increases. For the simulation examples in this paper, we focus on an exponential covariance structure, $∑(σ,κ)$ where $σ=(σ1,…,σp)$ and $κ>0$⁠. The element, $∑ij$⁠, in the ith row and jth column of this covariance is

Different parametric covariance structures give different forms for how the correlation between $yi$ and $yj$ decreases as $d(θi,θj)$⁠. For an exponential covariance, the parameter $κ$ controls how rapidly the contour can change over space, that is, the roughness or smoothness of the contour. For example, Figure 4 shows sample contours with the same $μ$ and $σ$ values as in Figure 3. However, the contours appear smoother because of the different covariance structure. The resolution of the data should also be considered in considering the roughness of the contour. If the data resolution is high, information about how smooth or rough the contours are can be assessed and modelled accurately with the $κ$ selection. If the resolution is low, though, only smoother contours will be possible, and $κ$ must be restricted accordingly. In other words, estimates of $κ$ from data may be limited by the data's resolution.

3.2 Fitting GSCMs

We now turn to building a GSCM given observed contours. We assume that the data are $N$ observed contours, $S=(S1,…,SN)$⁠, that enclose regions that are star-shaped polygons. (Note that each $S$ may be derived from binary grid as described in Section 2.) We also assume that the contours are generated from a common, but unknown, $C$ and $θ$⁠, that is, $C=C1=⋯=CN$ and $θ=θ1=⋯=θN$⁠. We first find a starting point, $C^$⁠, and angles, $θ^$⁠. Then, we estimate $μ$ and the parameters controlling $∑(·)$ based on the observed $y$ for the selected $C^$ and angles $θ^$⁠.

3.2.1 Fixing the starting point $C^$ and the set of angles $θ^$

The accuracy of the GSCM depends on the selection of $C^$ and $θ^$⁠. If $C^$ is selected outside the kernel of the polygons, the line set will not be able to intersect with all sections of the contour line. So some area within the true enclosed contour will be missed by the contour model. Similarly, if $θ^$ is not dense enough, some direction changes in the contour will not be modelled. The influence of the density of $θ^$ will be greater if the contour changes direction rapidly. The set of angles, $θ^$⁠, is selected to keep the mean difference in area between the observed contours' enclosed polygons and the star-shaped representations of these contours' enclosed polygons below some value.

Following Corollary 1, any set of distinct angles can be used to form a star-shaped representation of a contour. We recommend using evenly-spaced angles. Figure 5 illustrates how a star-shaped contour is approximated with a star-shaped representation with evenly-spaced $θ$⁠. Using evenly spaced angles ensures that the model represents all sections of the contours with the same level of precision. While some sections of the contours may change more rapidly than others and require more lines to represent well, exactly where these changes occur is hard to determine without having already modelled the contour. Therefore it is easier to select a density of lines to use consistently for all sections of the contour.

This sensitivity of GSCMs to the choice of $C^$ and $θ^$ motivates the development of a careful procedure to select them. Our procedure seeks to minimise the area that cannot be represented by the GSCM while balancing computational constraints. We fix a starting point, $C^$⁠, and set of angles, $θ^$⁠, that can be used to describe the observed contours accurately. We first describe how to find $C^$ conditional on $θ$ and $θ^$ conditional on $C$ separately. Then, we describe an iterative algorithm to fix both values together.

Finding$C^$conditional on$θ$⁠: The starting point $C^$ used in modelling is selected to minimise the difference in area between the observed contours' enclosed polygons and the star-shaped representations of the observed contours' enclosed polygons. Conditional on $θ$⁠, we define the set that differs between an observed polygon, $S_i$⁠, and a star-shaped representation of that contour with starting point, $C^$⁠, as

where the superscript $c$ denotes the complement of the set. The area contained within this set is denoted by $|A(C^,θ,Si)|$⁠. The lower this area, the more accurately the GSCM will be able to represent the enclosed polygons. So, selecting a $C^$ such that $|A(C^,θ,Si)|=0$ for all $i$ would be optimal. Assuming $θ$ is known, Theorem 1 guarantees that at least one point $C^$ exists where $|A(C^,θ,Si)|=0$ for all $i$⁠. If only one point satisfies this condition, then $C^=C$⁠. However, any location $C^$ for which $|A(C^,θ,Si)|$ is approximately zero will be able to represent the area within the contour well. As such, we select $C^$ numerically with

where $f$ is the mean function and $𝒟$ is a grid of points. The finer the grid of $𝒟$ the closer $|A(C^,θ,Si)|$ will be to zero. The number of points in $𝒟$ can be reduced by restricting assessment to the kernels of star-shaped polygons.

Theorem 2

Let$S_$be a set of$N$star-shaped polygons, let$C$be the true starting point, and let$𝒦^(S_)=𝒦(S_1)∩…∩𝒦(S_N)$denote the intersection of the kernels of all polygons. Then,$C∈𝒦^(S_)$⁠. (Proof in Appendix A.)

Therefore, we need only perform the optimisation in Equation (5) for $𝒟∈𝒦^(S_)$⁠. Algorithms for computing $𝒦^(S_)$ are given in Appendix  B.

Finding$θ^$conditional on$C$⁠: As motivated in Section 3.2.1, GSCMs can be fit with evenly spaced angles. So, finding the set of angles $θ^$ only requires finding $p^$⁠, the number of elements in $θ^$⁠.

Since $p^$ controls the dimensions of $μ$ and $∑$⁠, larger $p^$ requires more computation. We then identify the approximately lowest $p^$ that keeps the mean difference in area below some value. To allow comparisons of polygons of different sizes, we often express this allowable mean difference in area as a proportion, $δ$⁠, of the average area of the polygons. This constraint on the allowable differing area is then

where $f$ is the mean function and $|Si_|$ denotes the area of the polygon $S_i$⁠. We use Algorithm 1 to find an approximately minimal $p^$ that satisfies the constraint in Equation (6). To avoid selecting a larger $p^$ than necessary, $p(0)$ should generally be initialised such that Equation (6) is not satisfied. Smaller values of $a$ will make $p^$ more precise, but will require more computation than larger $a$⁠. Using lower values of $δ$ will generally result in higher $p^$ and lower differences in area. Section 5.3 uses simulation to explore different $δ$⁠.

Finding$C^$and$θ^$⁠: In practice, neither $C^$ nor $θ^$ will be known. So, to find both values, we iterate between setting $C^$ conditional on $θ^$ and $θ^$ conditional on $C^$⁠. Algorithm 2 describes this process. As in Algorithm 1, the initial value, $p(0)$⁠, should be selected to be low enough that Equation (6) is not satisfied. Otherwise, we may select a larger $p^$ than is needed. Smaller values of $a$ will result in more precise determination of $p^$⁠, but will require more computation.

3.2.2 Computing a posterior

Once $C^$ and $θ^$ are determined, model fitting is straightforward. For each observed contour, the observed $y$ values are computed given $C^$ and $θ^$⁠. Then from Equation (2), the corresponding likelihood for these $y$ is just that of a multivariate normal distribution:

Estimates of the mean vector $μ$ and the parameters controlling $∑(·)$ can be estimated using any standard method.

For demonstration purposes in the simulations and for modelling the sea ice edge we take a Bayesian approach. We assume an exponential covariance as in Equation (3) with parameters $σ$ and $κ$⁠. We use the following simple prior distributions for $μ$⁠, $σ$ and $κ$⁠. We assume a multivariate normal prior distribution for $μ$⁠:

The hyperparameters $μ0$ and $Λ0$ give the mean and covariance of the prior distribution. For the example exponential covariance defined in Equation (3) we assume uniform priors on $κ$ and $σ$⁠:

where the hyperparameters $βκ,0$ and $βσ,0$ are upper bounds. Samples from the posterior distributions of the parameter can be found via standard Markov chain Monte Carlo (MCMC).

3.2.3 Estimating gridded probabilities and credible intervals

Sampled contours can be used to estimate the probability of a given area being contained within a contour. Consider $K$ sample contours, $S^=(S^1,…S^K)$⁠. Each generated contour can be approximated by a binary grid, $G$⁠, of dimension $r×v$⁠. Let $gi,j$ indicate the grid box in the ith row and jth column of $G$⁠. Let $1gi,j,k$ indicate whether the majority of the area in grid box $gi,j$ is inside or outside contour $k$⁠. Most grid boxes will be entirely inside or outside the contour; however, grid boxes that intersect the generated contour will contain area both inside and outside the contour. Ideally, the grid selected should be fine enough to ensure that little area is contained within these transitional grid boxes. Averaging the binary grids produces an $r×v$ matrix, $P^$⁠, with elements $p^i,j=∑i=1K1gi,j,k/K$⁠, that indicate the probability of grid box $gi,j$ being contained within a contour. The last row of Figure 3 shows estimated gridded probabilities obtained from $K=100$ generated contours from the corresponding GSCMs.

Credible regions for the location of the contour can be computed from $P^$⁠. The $(1−α)$ credible region, $I1−α$⁠, is formed from a union of grid boxes that satisfy the condition

3.2.4 Rescaling data

For numerical convenience in fitting and generating contours, it is often desirable for all contours to be contained within the $[0,1]×[0,1]$ unit square. Observed data will typically need to be rescaled to be within these bounds. Data should be re-scaled such that generated contours do not extend outside the unit square. A good re-scaling also ensures that the contours that will be generated rarely, if ever, extend outside the unit square.

Therefore, we re-scale observed contours $S=(S1,…SN)$ to be within an $[ϵ,1−ϵ]×[ϵ,1−ϵ]$ square. This re-scaling provides a buffer region of width $ϵ$ on the outside of the unit square in which no contours have been observed. Therefore, if generated contours extend farther than the observed contours, they will typically go into this buffer region rather than outside the unit square. The higher the variability of contours, the larger the value of $ϵ$ needed to avoid generating contours that go beyond the unit square.

To transform a set of observed coordinates, $S=(S1,…,Sn)$⁠, to the square of dimension $[ϵ,1−ϵ]×[ϵ,1−ϵ]$⁠, let $min(Sx)$ and $max(Sx)$ denote the minimum and maximum observed $x$-coordinates from all spatial points in all contours in $S$⁠. Define $min(Sy)$ and $max(Sy)$ analogously for the $y$-coordinates. Let $δ=max(max(Sx)−min(Sx),max(Sy)−min(Sy)))$ and let $si,j=(si,jx,si,jy)$ denote the ith point in the jth.

Then for all $i$ and all $j$⁠, we shift and re-scale all coordinates for observed points $si,j$ with the transformation

Note that Equations (12) and (13) adjust the coordinates such that all lengths $y$'s are re-scaled by the same constant and so still follow a Gaussian distribution.

4 COVERAGE METRIC

To assess if our probabilistic contour model performs well in representing the ice edge, a metric is needed. A good model correctly identifies the region where the contour could plausibly be located. So we focus on the coverage of prediction intervals for star-shaped contours. With an accurate contour model, the variability of the generated contour would be correctly represented along all parts of the contour. In designing an appropriate metric, we leverage the star-shaped structure of the data. The general idea is to assess coverage for each line in a line set individually.

To make this idea precise, we define several quantities illustrated in Figure 6. As in Equation (11), let $I1−α$ be the $1−α$ credible region obtained from some contour model. Define some test line set $ℒ∗(C∗,θ∗)$ with $M$ evenly spaced lines. Define $I1−α,k$ as the line segment formed from the intersection of the $I1−α$ credible region and the line $ℓk∈ℒ∗(C∗,θ∗)$⁠. We refer to $I1−α,k$ as a test line. Also, define $Ri,k$ as the intersection of some observed contour $Si$ and line $ℓk$⁠. Note that $Ri,k$ will always be a single point when the polygon $S_$ is exactly star-shaped and $C∗∈𝒦(Si_)$⁠.

Let $Wi,k=1[Ri,k∈I1−α,k]$ indicate whether the intersection points of the observed contour and the line $ℓk∈ℒ∗(C∗,θ∗)$ are contained within the intersection of the credible region and line $ℓk$⁠. Then, for credible intervals with perfect coverage, for any $i,k$⁠,

In other words, the part of the true contour that intersects line $ℓk$ is contained within the credible contour region along line $ℓk$ with probability $1−α$⁠.

We consider coverage behaviour for a set of $N$ observed contours, $S=(S1,…,SN)$ that enclose star-shaped polygons assumed to have been generated from the same model. Each $Si$ is assumed to be independent of $Sj$ for all $i,j$⁠.

For a credible interval with perfect coverage for any $i$⁠, $k$⁠,

for sufficiently large sample size, $N$⁠. This equation is used to assess coverage in practice. Equation (15) not holding for any $k$ indicates the variability of the contour on test line $ℓk$ is not correctly represented by the credible region.

With this set-up, $Wi,k$ and $Wj,k$ are independent for all $i,j$ conditional on the parameters of the GSCM. So for a given line $ℓk$⁠,

Since the distributional behaviour of the quantity $∑i=1NWi,k$ is known, the expected variability of the mean coverage given a particular sample size is known. This information can be accounted for in a simulation study or cross-validation experiment.

We also consider how coverage on one test line in $ℒ∗(C∗,θ∗)$ relates to coverage on another test line. The location of points on a contour are generally correlated. So, $Wi,k$ and $Wi,ℓ$ are correlated. Typically $Wi,k$ and $Wi,ℓ$ will be more correlated the smaller the angle distance from $θk$ and $θℓ$⁠. What this means is that while

is fixed and known, the quantity $∑k=1MWi,k$ is not a good metric to assess to coverage. The distribution of $∑k=1MWi,k$ depends on the correlation structure of the points on the contour. For a contour with high correlation among points, the quantity $∑k=1MWi,k$ will be substantially affected by what contour happens to be sampled. For intuition, consider a contour with high correlation among all $i,j$⁠. In this case, the contour is likely to be entirely within the credible interval or entirely outside of it. So, $∑k=1MWi,k$ will be either 0 or M.

These relationships among coverage for $M$ and $N$ show that sample sizes need to be considered in terms of the number of contours observed. The metric $∑i=1NWi,k$ should be considered for all elements in $θ$⁠. The number of test lines, $M$⁠, should be set such that accuracy is assessed with detail appropriate for the application. Since the true exact value of $p$⁠, is unknown, we cannot simply set $M=p$⁠. However, based on the observed data, we should have a general idea of $p$⁠. So, we set $M>>p$⁠, to ensure that $M>p$⁠. We also evenly space $θ∗$ with $θ1∗=(2π/M)/2=π/M$⁠. By assessing coverage on a substantially greater number of test lines than the true number of lines, we ensure that coverage is at least assessed near every true line.

To carry out this assessment, a fixed starting point $C∗$ and $θ∗$ should be selected. In simulation studies, the true value of $C$ will be known and we can let $C∗=C$⁠. For assessment of real data such as in a cross-validation study, a starting point for $C∗$ will be unknown and must be determined. We recommend using a $C∗$ that minimises the difference in area between the observed contours' enclosed polygons and the star-shaped representations of the observed contours' enclosed polygons as in Section 3.2.1, that is, let $C∗=C^$⁠.

This assessment approach differs from how contours have been assessed in the context of level exceedances. There, a credible region or confidence region has often been defined as the region that covers the true contour in its entirety $(1−α)$-proportion of the time (e.g., Bolin & Lindgren, 2015; French, 2014). With credible (confidence) regions constructed to satisfy this definition, coverage can be assessed by determining what proportion of the time the true contour is fully contained within the region. We opt not to use this metric since our goal is to develop a method to generate ice edge contours directly. Correlation along the contour makes assessing the probability of capturing the entire contour difficult. Our metric reflects that we are most concerned with getting the right variability in all parts of the contour. We are less concerned with identifying a larger area that contains that entirety of the contour with high probability. Our intervals are therefore narrower than would be required for these global intervals.

5 SIMULATION STUDIES

5.1 Simulation details

Before evaluating the ice edge data, we consider how the star-shaped model performs in inferring distributions of simulated data. We consider performance with varying numbers of observations, different constraints for the allowable mean difference in area (⁠$δ$ as defined in Section 3.2.1), and varied GSCM parameters.

In many of our simulations, we focus on a particular GSCM with $p=50$ that we will refer to as Shape A. The correlation structure of $y$ follows the exponential form given in Equation (3). The vector of mean distances, $μ$⁠, and variance parameter vector, $σ$⁠, change gradually. The exact values of $μ$ and $σ$ can be found in Appendix C. Unless otherwise noted, $κ$ is set to 2. Example generated contours and gridded probability estimates for Shape A are given in the left panel of Figure 3.

The values $C^$ and $θ^$ are found as described in Section 3.2. The parameter values are fit from observations using MCMC. Chains are run for 50,000 iterations with the first 15,000 iterations discarded as burn-in. The prior parameters are $μ0=(0.2,…,0.2)$⁠, $βκ0=(8,…,8)$⁠, $βσ0=(0.15,…,0.15)$⁠, and $Λ0=.05Ip^$ where $Ip^$ is a diagonal matrix of dimension $p^$ by $p^$⁠. The value $p^$ refers to the number of angles in $θ^$⁠.

For all simulations, we use 40 evaluation runs. On each evaluation run, we estimate the GSCM parameters using $N$ contours generated from the true GSCM as training data. From the resulting fitted GSCM, we generate 100 contours and find credible intervals as described in Section 3.2.3. To evaluate coverage, we compare the credible interval with a single ‘true’ contour drawn from the true GSCM. We record the coverage for a set of $M=100$ evenly spaced test lines with $θ1∗=π/M$ as described in Section 4. We report the mean coverage over the 40 evaluations and the SD across the $M=100$ test lines.

5.2 Varying number of observations, $N$

Our first simulation varies $N$⁠, the number of simulated ‘observed’ contours used to fit the Shape A GSCM. We consider coverage performance for 10, 20 and 50 simulated observed contours with $δ=0.02$⁠. Table 1 displays the results of this simulation.

We find that coverage improves for $N=20$ compared to $N=10$⁠. We find slightly worse performance for $N=50$ than $N=20$⁠, although the performances are not substantially different given that we had only 40 evaluation runs. These results indicate that obtaining some minimum sample size is important for coverage performance. For small sample sizes, on the order of $N=10$⁠, the data alone may not be enough to produce accurate coverage, particularly for the 80$%$ credible interval. However, data sets of this size can potentially still be modelled correctly if informative priors can supplement the observations.

5.3 Varying allowable difference in area, $δ$

In these simulations, we evaluate how coverage accuracy is affected by the parameter $δ$ introduced in Section 3.2.1. This parameter controls the allowable differing area when setting the number of lines used in fitting, $p^$⁠, and how accurate $C^$ must be. We evaluate coverage for Shape A with $δ$ set to $0.03,0.02$ and 0.01. These $δ$ selections set the allowable mean difference in area to 3$%$⁠, $2%$ and $1%$ of the mean area contained within the observed contours. We also consider how correlation in $y$ affects the need for different $p$ by evaluating each $δ$ for three different $κ$ values: 1, 2 and 4. Table 2 displays the mean coverage across test lines for three $α$-levels along with the mean $p^$ found. On each evaluation the number of sampled contours is set to $N=20$⁠.

We find that the mean coverage accuracy is only modestly affected by the value of $δ$⁠. These results support the idea that using a lower $δ$ in many cases will reduce computation while not reducing model performance. We also find that, for a given $δ$⁠, an increase in $κ$ corresponds to a decrease in $p^$⁠. In other words, for a contour with higher correlation among $y$⁠, a smaller set of lines can adequately represent the contour distribution.

5.4 Varying GSCM parameters

We also evaluate contour models that have mean values, $μ$⁠, that vary more slowly and more quickly than in Shape A. These GSCMs are denoted as Shape B and Shape C, respectively. Figure 3 shows sample contours and probability distributions for these shapes. Both models are defined to have $p=50$ and $κ=2.$ Exact values for $μ$ and $σ$ for Shape B and Shape C can be found in Appendix C. We report coverage results in Table 3 for three $α$-levels using $N=20$ simulated observed samples and $δ=0.02$⁠.

We find reasonably accurate coverage performance for all three shapes. The method performs slightly worse for Shape C than for Shape A and Shape B. This performance difference is likely due to the difficulty of getting the pointed sections of Shape A in the correct location if $p^$ is even slightly underestimated. This result indicates that for contours that look like Shape A, lower $δ$ values may be appropriate. We note that the simulations in this section are relatively symmetric. Results might be more variable moving around the contours if there was less symmetry. In such cases a lower $δ$ would be appropriate to account for this heterogeneity. Overall, the good performance across parameter settings indicates that a range of contours can be well approximated by a GSCM.

6 MODELLING CONTOURS ENCLOSING APPROXIMATELY STAR-SHAPED POLYGONS

The ice edge contours enclose polygons that are approximately, but not exactly, star-shaped. This section describes how to assess whether the GSCM is appropriate given the observed data, and how the fitting procedure is altered if the observed contours enclose polygons that are not exactly star-shaped.

6.1 Assessing appropriateness of GSCM

Two main assumptions must be met to apply the GSCM: the polygons enclosed by the contours must be approximately star-shaped and all contours should have at least one common point. The latter assumption is needed to define a starting point and can be trivially assessed. The former assumption can be assessed using metrics that describe how close an observed contour is to enclosing a polygon that is star-shaped. These metrics focus on the difference in the area between the polygon enclosed by the true contour and the polygon enclosed by the star-shaped representation of the contour. If these differences are small for a set of observed contours, $S$⁠, then the GSCM can be applied.

We relax Definition 8 to make precise how to approximate an arbitrary polygon with a star-shaped representation. Two main differences between star-shaped polygons and arbitrary polygons are addressed in these new definitions: (1) an arbitrary polygon may not have a kernel and (2) the contour line enclosing an arbitrary polygon may intersect with some of the lines in the line set multiple times. The new star-shaped representation definitions are:

Definition 10

Underestimated star-shaped representation, $V˜u(C,θ,S)$⁠: Let $S_$ be a polygon described by ordered spatial points $S=(s1,…,sn)$⁠, let $C∈S_$ be a starting point, let $θ$ be an arbitrary set of $p$ unique angles, and let $y=(y1,…,yp)$ be a set of distances from $C$ to the closest intersection point of the contour line $S‾$ and each line $ℓi$ in the line set $ℒ(C,θ)$⁠. Then, the star-shaped representation of the contour, $V(C,θ,y)$⁠, is the underestimated star-shaped representation, $V˜u(C,θ,S)$⁠.

Definition 11

Overestimated star-shaped representation, $V˜o(C,θ,S)$⁠: Let $S_$ be a polygon described by ordered spatial points $S=(s1,…,sn)$⁠, let $C∈S_$ be a starting point, let $θ$ be an arbitrary set of $p$ unique angles, and let $y=(y1,…,yp)$ be the set of distances from $C$ to the farthest intersection point of the contour line $S‾$ and each line $ℓi$ in the line set $ℒ(C,θ)$⁠. Then, the star-shaped representation of the contour, $V(C,θ,y)$⁠, is the overestimated star-shaped representation, $V˜o(C,θ,S)$⁠.

The names of these representations highlight that these polygons generally under- or overestimate the area contained within the true polygon. For contours that enclose star-shaped polygons, if $C∈𝒦(S_)$⁠, only one intersection is found between the contour line and all lines in the line set. So, $V˜o(C,θ,S)=V˜u(C,θ,S)=V˜(C,θ,S)$⁠. For notational convenience, let $V˜u=V˜u(C,θ,S)$ and $V˜o=V˜o(C,θ,S)$⁠. Then the sets that differ between the true contour and the under- and overestimated star-shaped representations are

where the superscript $c$ denotes the complement of the set. Let $|Au(C,θ,S)|$ and $|Ao(C,θ,S)|$ be the area contained within these sets. Figure 7 illustrates the quantities described in this section for four contours. The difference in area is zero only if the polygons are star-shaped. More precisely:

Theorem 3

For any polygon$S_$that is not star-shaped$|Au(C,θ,S)|>0$and$|Ao(C,θ,S)|>0$for any$C$and$θ$⁠. (Proof in Appendix A.)

For easier comparison of these differences in area across sets of polygons of different sizes, these areas can be expressed as percentages of the mean total area of the polygons. Table 4 reports the difference in area for the contours in Figure 7, and their star-shaped representations as a percentage of the total area of the polygon.

6.2 Fitting GSCMs to approximately star-shaped polygons

The approach to fitting in Section 3.2 needs to be altered slightly for contours that enclose regions that are only approximately star-shaped contours. The values of $y$ need to be computed using the under- or overestimated star-shaped approximation as given in Definitions 10 and 11. Whether to use the under- or overestimated star-shaped representation depends on the application. In some cases, asymmetric risks may motivate selecting a model that generally over- or underestimates the area within the polygon. Otherwise, both $|Au(C,θ,S)|$ and $|Ao(C,θ,S)|$ can be computed for the set of observed contours $S$ and whichever representation results in less difference in area can be selected. Fitting the posterior proceeds as in Section 3.2.2. Finding probabilities and credible intervals proceeds as in Section 3.2.3.

We find $C^$ and $θ^$ using nearly the same algorithm as in Section 3.2.1, except that we update the star-shaped representation in Equations (5) and (6) to be the under- or overestimated star-shaped representation. Specifically, we replace Equation (5) with

for the underestimated representation and

for the overestimated representation. The values $C^u$ and $C^o$ denote the estimated starting points for the two approximations, respectively. As before, $f$ is the mean function. Like in Section 3.2.1, we evaluate these functions of the area difference at a grid of possible locations, $𝒟$⁠. Unlike in Section 3.2.1 we cannot use the estimated intersection kernel to reduce the size of $𝒟$⁠, since the polygons are not exactly star-shaped. However, we note that the starting point must be common to all observed polygons, so we constrain $𝒟$ to only cover areas in the intersection of all observed polygons, $S_1∩…∩S_N$⁠. Equation (6) is also slightly changed so that

for the underestimated representation and

for the overestimated representation. The value $f$ still represents the mean function and $δ$ is still a proportion.

6.3 Coverage metric for approximately star-shaped polygons

This metric introduced in Section 6.1 is essentially the same when applied to approximately star-shaped polygons. However, $Ri,k$⁠, the intersection of an observed contour $Si$ and line $ℓk$⁠, may now contain multiple points since polygon $S_i$ is only approximately star-shaped. A single point of intersection will still be more common when polygons are approximately star-shaped. Aside from this change, the definition of $Wi,k=1[Ri,k∈I1−α,k]$ remains the same when $Ri,k$ is composed of multiple points. So, all subsequent definitions and properties in Section 6.1 hold as well.

6.4 Simulation study of contours enclosing approximately star-shaped polygons

In these simulations, we assess GSCM performance for contours that enclose polygons that are only approximately star-shaped. We simulate contours that vary systematically in how much the polygons they enclose differ from being star-shaped. Figure 8 shows examples of the types of contours we will evaluate.

To obtain these contours, we first generate polygons that are star-shaped from a GSCM. We then append sections to these polygons that cause the polygons to no longer be star-shaped. The appended sections loop back around the outside of the initial polygon over some number of lines in the initial line set. The number of initial lines looped back over are selected randomly from some uniform distribution. Appended sections that loop around a larger number of lines are longer than appended sections that loop over a smaller number of lines. Longer appended sections result in more area that cannot be described with a star-shaped representation than shorter appended sections. How close the appended section is to the initial star-shaped polygon and the width of the appended section are selected randomly from uniform distributions.

We consider 40 evaluation runs for three different uniform distributions for the number of initial lines that the appended section loops back over. The initial GSCM is set to be Shape A with $κ=2$ and $p=50$⁠. In each evaluation run, the GSCM is fit to $N=20$ simulated contours. Rather than fixing $δ$ to determine $p^$ and $C^$ in these simulations, we set $p^=p=50$⁠. We then estimate $C^$ conditional on $p^$⁠. This choice simplifies the interpretation of our results. (Larger $δ$'s would be needed as the area differing from a star-shaped polygon increases. So, it would be difficult to distinguish whether performance differences were due to changes in $p^$ or changes in how close to star-shaped the polygons are.) We initially run these simulations using a random location for the appended section and then repeat these simulations with a fixed location for the appended section. Results are reported in Table 5.

With a random location for the appended section, we find that applying GSCMs still results in reasonable coverage for the 90$%$ and $95%$ credible intervals. For the $80%$ credible interval, performance is moderately degraded. Interestingly, average performance does not seem to be correlated with how large the appended section to the contour is. When an appended section is added to a fixed location, we find that the mean coverage across the test lines is relatively accurate; however, the SD is quite high, suggesting that coverage is actually poor in some parts of the contours. Figure 9 illustrates this variability in performance for the case with the number of initial lines looped over distributed Uniform(4, 5). We plot the proportion of evaluation runs covered by the $90%$ credible intervals for each test line individually.

The location of the fixed appended section is under-covered. In contrast, no obvious patterns are seen when the location of the appended section is random. These results indicate that contours that enclose polygons that modestly differ from star-shaped contours can be modelled with GSCMs. However, if areas differing from the star-shaped representation occur in the same location repeatedly, additional modelling of these areas may be needed to avoid systematic errors in coverage.

7 APPLYING GSCMS TO THE ARCTIC SEA ICE EDGE CONTOUR

We now illustrate how GSCMs can be used to model the Arctic sea ice edge. We focus on the ice edge in September, the month when the ice-covered area in the Arctic is at its annual minimum. September holds particular interest for maritime planning, since vessel traffic is typically highest when there is the least sea ice.

We focus on modelling the sea ice edge over a short set of recent years. These types of models can be used to make statements about how probable sea ice edges will be in the near future and/or to describe plausible sea ice edges that align with the oceanic and meteorological conditions in this time period. The GSCM can be applied to generate plausible sea ice edge contours and corresponding credible intervals for the ice edge contour. Note that with climate change, the Arctic area covered by sea ice (Comiso et al., 2008; Stroeve et al., 2012) has reduced on average over time. So, observed ice edge contours from decades past are not similar to ice edges observed recently. However, on short time scales, the change in average meteorological and oceanic conditions is small. So, assuming that each observed September sea ice edge is an independent draw from a stationary distribution is appropriate. Additionally, the correlation of sea ice from 1 month to the next decays rapidly, so observations from one year apart can be treated as independent.

The data used in fitting a sea ice edge model is a monthly average observational product produced by the National Aeronautics and Space Administration satellites Nimbus-7 SMMR and DMSP SSM/I-SSMIS and downloaded from the National Snow and Ice Data Center (Comiso, 2017). The data are composed of gridded fields reporting the sea ice concentration, or percent of ice-covered area in each grid box. Following a convention used by sea ice researchers, grid boxes with at least 15$%$ concentration are treated as containing sea ice and grid boxes with less than 15$%$ concentration are treated as open water. This thresholding is needed because satellite estimates for very low concentrations are not considered reliable. The transition from complete sea ice cover to open water occurs over a narrow spatial range, so the ice edge is not particularly sensitive to the exact threshold selected. Each month's gridded field is converted to an ordered sequence of points as described in Section 2. We refer to the ordered sequences of points associated with each month of data as the observed ice edge in that month.

To assess the sea ice edges generated by the GSCM, we perform a leave-one-out cross-validation experiment on the September sea ice edge contour in a region in the central Arctic for 10 recent years. For each year $j$ from 2008 to 2017, we fit the GSCM using the contours observed in the other 9 years. We then try to ‘predict’ the distribution of possible ice edges that would have been plausible in year $j$⁠. Data have been re-scaled as described in Section 3.2.4 with $ϵ=0.1$⁠. We exclude from fitting a section of the ice edge contours that always borders land, since there is no variability in the associated $y$ values. We use the overestimated approximation, since for some maritime planning applications, predicting too much sea ice is less dangerous than predicting too little sea ice.

As shown in Section 5.2, 9 years of data may not be sufficient for fitting a GSCM well. However, we have considerable information that we can incorporate into the prior to augment the observations. In particular, the region we are focused often follows land boundaries, restricting where sea ice in the region could go. So, the lines in the line set $ℒ=(ℓ1,…,ℓp)$ are also bounded. We leverage their fixed maximal lengths in setting the priors. Specifically, we let $μ0,i=0.5ℓi$ and $β0,σi=(ℓi/2)/Φ−1(.995)$⁠. The latter corresponds to the standard deviation of a normal distribution with $99%$ of its mass falling in the interval $(0,‖ℓi‖)$ where $‖ℓi‖$ is the length of $ℓi$⁠. We also set $Λ0$ be a diagonal matrix with $0.52$ on the diagonal and $βκ0=10$⁠. We use MCMC for fitting with 50,000 iterations, of which 25,000 are omitted as burn-in.

We evaluate coverage as in Section 4 using the approximately optimal starting point selected using the method in Section 6.2 with $δ=0.03$ and $p=100$⁠. We exclude from evaluation test lines that always border land and have no variability. This leaves 68 test lines. We plot in Figure 10 the 80$%$ credible interval estimated for the 2017 sea ice edge fit with data from the other nine years. As in Section 4 lines extending outward from the estimated starting point are coloured blue if they do not cover the left-out contour and black otherwise. Table 6 reports the mean coverage of the credible intervals averaged over the 10 years and 68 testing lines. Three nominal $α$-levels are considered. The observed and nominal coverage values are similar, suggesting that the GSCM has appropriate coverage. In other words, GSCMs perform well for capturing the variability of the sea ice edge contour.

8 DISCUSSION

We have introduced the GSCM for modelling contours that enclose polygons that are star-shaped or approximately star-shaped and we have applied it to ice edge contours. Analysis of September Arctic sea ice showed how GSCMs can be applied to model the ice-covered area in the Arctic. Simulation studies also illustrated how GSCMs provide accurate coverage in different scenarios. We conclude this paper with a discussion of how GSCMs relate to other contour models, other methods' appropriateness in modelling the sea ice edge, and directions for future research.

8.1 Other approaches

A large body of research addresses contours in other contexts. GSCMs are applied when multiple contour boundaries are directly observed and the distribution of possible contours is the primary object of inference. Other methods for contours and boundaries may be appropriate for other applications.

Much of the existing research on boundaries relates to level exceedances, also called excursions. The classic level exceedance problem refers to inferring the contour enclosing regions where some latent process exceeds a certain level $u$⁠. Inference is based on measurements taken at various spatial points on a random spatial field. Early work by Polfeldt (1999) considers how to make statements about the accuracy of contour maps in this context. Lindgren and Rychlik (1995) first define contour uncertainty regions using unions of crossing intervals, or line sections where transitions from below and above $u$ occur. More recently, Bolin and Lindgren (2015) introduce a method for inferring exceedance levels with irregularly-spatial measurements when the latent spatial field is Gaussian. Their approach provides a way to make global statements about the uncertainty of the full contour. Bolin and Lindgren (2017) then extend this method to estimate the uncertainty of multiple contours to produce contour maps with appropriate uncertainty estimates. Both methods leverage Integrated Nested Laplace Approximations for efficient computation and can be implemented with the excursions R package (Bolin & Lindgren, 2018).

French (2014) provides an alternate simulation-based method for making global statements about the location of the contour. Methods for identifying the exceedance region are also explored from both Bayesian and frequentist perspectives (French & Hoeting, 2016; French & Sain, 2013). Level exceedance methods and GSCMs both focus on inferring contours and their uncertainty. However, in the former, boundaries and their uncertainty are inferred with measurements of a continuous process made at spatially referenced points while in the latter distributions of plausible contours are inferred from direct observations of the contour boundaries themselves.

Wombling methods also focus on spatial boundaries. First considered by Womble (1951), these methods typically apply bilinear interpolation to spatially-referenced data points. The gradients of the interpolated functions are used to infer boundaries. Jacquez et al. (2000) summarise early research on primarily deterministic methods for identifying these boundaries. Principled Bayesian statistical methods have since been developed (Banerjee & Gelfand, 2006; Gelfand & Banerjee, 2015) and recent research has introduced Wombling methods for areal data (Li et al., 2015; Lu & Carlin, 2005) and point processes (Liang et al., 2009). Wombling boundaries are inferred from spatially referenced or areal data. The proposed GSCM differs from Wombling techniques, since it is targeted to be applied to repeated directly observed contour boundaries.

Statistical shape analysis (e.g., Dryden & Mardia, 2016; Srivastava & Klassen, 2016) describes features and variation around boundaries. Shapes typically have consistent and definable features. In these types of applications, location and rotational effects are often ignored in describing the distribution around the shape. As an example, shape analysis is often applied to biological imaging research. Deformable templates were developed to describe distributions around shapes with definable features such as the parts of a hand (Amit et al., 1991; Grenander & Keenan, 1993; Grenander & Miller, 1998). A key difference between our application and those more typical of shape analysis is that the sea ice edge's physical location if of primary importance. As such, establishing correspondence (Procrustes analysis) is not appropriate. We are interested in the specific location of the ice edge on Earth, not the general shapes that form ice edges. As such, translating, rotating, or scaling observed points to align one ice edge to another does not make sense. Furthermore, no discernible features are present in the ice edge to align.

Many image analysis methods have been developed to segment images or identify edges. Mathematical morphology (Haralick et al., 1987; Lee et al., 1987), watershed segmentation (e.g., Gauch, 1999), Bayesian models (Li & Ghosal, 2017) and more recently deep learning have all been applied to identify or sharpen the uncertainty of a single observed boundary in an image. The goals of these methods again differ from the goal of GSCM to define variability over multiple observations of boundaries.

Another alternative to GSCMs would be to model directly whether the points on a lattice are inside or outside a contour boundary. However, methods for modelling binary data on a lattice such as the autologistic (Besag, 1974), centred autologistic (Caragea & Kaiser, 2009) and the spatial generalised linear mixed model (Besag et al., 1991; Diggle et al., 1998; Hughes & Haran, 2013) are not structured to guarantee that all the grid boxes inside the contour form in a contiguous section. Hence, these methods are not designed directly for modelling contour boundaries.

Extensive research in the remote sensing community has focused on how best to translate the sea ice backscatter from Synthetic Aperture Radar into sea ice classifications. Approaches include thresholding of the backscatter coefficient for particular ice types or thickness, traditional statistical approaches and machine learning methods (see Zakhvatkina et al., 2019 for a thorough review). This research estimates the location of the ice edge given the corresponding satellite data for that time point. In other words, this research, like much of the statistical research on contours, focuses on finding the best estimate and uncertainty of individual contours at individual time points given the data at those time points. Instead, the GSCM allows one to infer the distribution of the ice edge given samples of ice edges observed at different times. This provides a distribution of plausible ice edges and information about the expected ice edge location and its variability.

Other sea ice research has examined the variability of the sea ice edge over time by parameterising the ice edge as a curve and considering deviations from the mean parameterised curve (Divine & Dick, 2006; Shapiro et al., 2003). This approach is reasonable for characterising past variability, but may not be optimal for generating a distribution of contours. Particularly if variability is high, a generated boundary may not enclose a single polygon or may intersect itself often.

8.2 Fractal contours and GSCMs

We have treated contours as connected sequences of points, but many contours have fractal-like properties. We now discuss how a fractal contour could be converted to a connected sequence of points for modelling with a GSCM. A true fractal contour, represented by a set $F$⁠, could be approximately represented by a smaller set of points $S$⁠. In a 2-D Euclidean space, $ℝ2$⁠, consider a countable or finite set ${U}$ of circles of radius $δ$⁠. We say ${U}$ covers $F$ if $F⊂∪i=1∞Ui$ and we refer to the set ${U}$ as a $δ$-cover of $F$⁠. The value $Nδ(F)$ is the smallest number of circles of radius $δ$ that could be used to cover $F$ (Falconer, 2004, pp. 27-28). Let ${U∗}$ denote one covering that contains $Nδ(F)$ circles. Since the contour $F$ is assumed to be finite, the number of circles in ${U∗}$ will also be finite. In Figure 11, we plot a $δ$-cover over a visualisation of a fractal contour $F$ with several finite self-similar layers. The $δ$-cover plotted is for visualisation and may not contain exactly $Nδ(F)$ circles.

We can define the elements in the sequence of points forming the contour, $S$⁠, to be the starting points, ${M∗}$⁠, of the circles in ${U∗}$⁠. The starting points should be arranged in the order in which they would be touched if one were to trace over the fractal contour line, $F.$ Since the contour is a closed loop, where and in what direction to start tracing the fractal contour, $F$⁠, only affects the indexing of the starting points and not the contour line formed by connecting these points. With this procedure the distance from any point $F$ to a point in $S$ is no more than $δ$⁠. In general there are multiple $δ$-covers that contains $Nδ(F)$ circles; therefore, a criterion needs to be specified as to which set ${U∗}$ is used to define $S$⁠. For example, the set ${U∗}$ could be selected to be the $δ$-cover with $Nδ(F)$ circles that has the circle centre closest to the highest x- and highest y-coordinate in the domain, that is, the top right corner of the domain.

Once $S$ has been obtained, GSCMs can be used. In particular, the under- or overestimated star-shaped representation to any fractal contour can be found for $S$ as described in Definitions 10 and 11. A star-shaped representation can represent the area contained within some fractal contours fairly well. For example, the visualised $F$ and the overestimated star-shaped representation for $S$ obtained using the visualised $δ$-cover are similar. The differing area represents only 7.4% of the total area of the visualised fractal contour when $p=200$ lines are used in the line set.

8.3 Limitations and extensions

A limitation of GSCMs is that they can be applied only to contours that enclose polygons that are star-shaped or approximately star-shaped. Using multiple starting points would be one way to relax this assumption. From each starting point a different star-shaped contour could be generated. These individual contours could be combined to produce an overall contour that encloses a polygon that is not star-shaped as in the ‘divide and conquer strategy’ of Li and Ghosal (2017). Another promising direction for relaxing the star-shaped assumption is to develop a method for appending sections to an initial contour that encloses star-shaped polygons. Extending the way contours were generated in the simulation in Section 6.4 might provide an initial approach to do this.

Another limitation is that the GSCM assumes that all generated line lengths $y$ are drawn from a Gaussian distribution. While GSCMs can cover a range of shapes, in some applications the observed $y$ may be skewed, bounded above, or otherwise differ from being distributed approximately normally. Therefore, exploring alternatives to GSCMs that allow for more flexible distributions of $y$ would be valuable.

Our development of the GSCM assumes a single data source defining the contours. However, in many applications multiple data sources may provide information about the ice edge contour or sections of the ice edge contour. For example, the satellite data we used in testing GSCMs is not the only source of sea ice data. Multiple satellites operating over different locations and at different spatial and temporal scales track sea ice behaviour. Determining how to combine these data sources into a probabilistic contour model could well increase applications to which the GSCM could be applied and improve the accuracy of existing models.

All models introduced in this paper assume constant fractal properties, that is, the smoothness (roughness) of the contour is the same for the whole contour. However, smoothness can vary spatially. Extending the GSCM to account for this type of non-stationarity would broaden the range of applications. Research exploring multi-fractal star-shaped objects with locally varying parameters (Emery, & Alegría, 2020) could provide a path to making this extension.

Overall, the GSCM provides a promising avenue for modelling data composed of multiple observed contours. Representing contours with sequences of points, which are of lower dimension than spatial fields, could enable more detailed and efficient modelling of contour boundaries.

ACKNOWLEDGEMENTS

We thank the Editor, Associate Editor, and two anonymous referees for helpful comments. We thank Cecilia Bitz, Donald Percival, and Daniel Pollack for helpful discussions. Results were generated using the ContouR R package (Director & Raftery, 2020) and scripts accessible at https://github.com/hdirector/contourPaperScripts. Contributions by Hannah M. Director and Adrian E. Raftery were supported by NOAA's Climate Program Office, Climate Variability and Predictability Program through grant NA15OAR4310161. Hannah M. Director's contribution to this work was also based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1256082. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

DATA AVAILABILITY STATEMENT

The Bootstrap Sea Ice Concentration Data that support the findings of this study are openly available in the NASA National Snow and Ice Data Center Distributed Active Archive Center at doi.org/10.5067/7Q8HCCWS4I0R (Comiso, 2017).

REFERENCES

APPENDIX A PROOFS

A.1 Proof of Theorem 1

Let $θ=(θ1,…,θp)$ with $θi<θi+1$ and $θi∈(0,2π)$ for all $i$⁠. For a star-shaped polygon $S_$ there exists $θ$ and $y$ such that $V˜(C,θ,y)=S$ for any $C∈𝒦(S_)$⁠.

Proof

Consider the point $si∈S$⁠. Since $C∈𝒦(S_)$⁠, there exists a line segment, $Csi‾$⁠, from $C$ to $si$ that is entirely contained within the polygon $S_$⁠. Let $θi=atan2(si,y−Cy,si,x−Cx)$ where atan2(⁠$b$⁠, $a$⁠) is the two-quadrant arctangent representing the angle between the positive $x$-axis and the line segment from the origin to point $(b,a)$⁠. Then the corresponding $ℓi$ in the line set covers the line segment $Csi‾$⁠. By construction, the line $ℓi$ intersects $si$⁠. The line $ℓi$ also cannot intersect any other points on $S‾$⁠, since the existence of such points would violate the assumption that $S$ is star-shaped. Define $yi=(si,y−Cy)2+(si,x−Cx)2$⁠, then $si=v^i$⁠. Repeat this construction of $θi$⁠, $ℓi$ and $yi$ for all $si$⁠. Then $si=v^i$ for all $i$⁠, hence $V˜(C,θ,y)=S$⁠.

A.1.1 Proof of Corollary 1

Let $ℓθ$ denote the line that extends infinitely outward from $C$ at angle $θ∈(0,2π)$ and that intersects $S‾$⁠. For any $θ∈(0,2π)$⁠, the line $ℓθ$ is distinct, that is, $ℓθ≠ℓθ′$ for any $θ,θ′$ such that $θ≠θ′$⁠.

Proof

Let $θi<θ<θi+1$⁠. The line $ℓθ$ intersects $S‾$ in the line segment $sisi+1‾$⁠, but not at either $si,si+1$⁠. Since this fact holds for any $θi,θi+1$⁠, each line $ℓθ$ must be distinct.

A.2 Proof of Theorem 2

Let $S_$ be a set of $N$ star-shaped polygons, let $C$ be the true starting point, and let $𝒦^(S_)=𝒦(S_1)∩…∩𝒦(S_N)$ denote the intersection of the kernels of all polygons. Then, $C∈𝒦^(S_)$⁠.

Proof

The starting point $C$ will be in every $𝒦(S_i)$ by the definition of a kernel. So, any point p that is in $𝒦(S_i)$ but not $𝒦(S_j)$ for any $i$⁠, $j$ is not $C$⁠. Hence, $C$ must be in $𝒦^(S_)$⁠.

A.3 Proof of Theorem 3

For any polygon $S_$ that is not star-shaped $|Au(C,θ,S)|>0$ and $|Ao(C,θ,S)|>0$ for any $C$ and $θ$⁠.

Proof

We show that $|Au(C,θ,S)|>0$⁠. The proof for $|Ao(C,θ,S)|>0$ is analogous. The quantity $|Au(C,θ,S)|=0$ only if $(S_c∩V˜u_)=∅$ and $(S_∩V˜u_c)=∅$⁠. These sets are both empty only if $V˜u=S$⁠. Hence, we need show only that no $θ$⁠, $C$⁠, $y$ combination exists that allows $V˜u=S$⁠. Polygon $S_$ is not star-shaped. So, there exists at least one point $s∗∈S‾$ such that for any point $C∈S_$⁠, the line $Cs∗‾$ goes outside polygon $S_$⁠. Since the line $Cs∗‾$ exits $S_$ before reaching $s∗$⁠, there is at least one additional intersection point between line $Cs∗‾$ and contour line $S‾$⁠. Let $s∗∗$ denote this intersection. The lines $Cs∗‾$ and $Cs∗∗‾$ are at the same angle, denoted $θ∗$⁠. For either $s∗$ or $s∗∗$ to be $V˜u$⁠, we must put a a line $ℓ∗$ in the line set that extends at angle $θ∗=atan2(sy∗−Cy,sx∗−Cx)=atan2(sy∗∗−Cy,sx∗∗−Cx)$ where atan2(⁠$b$⁠, $a$⁠) is the two-quadrant arctangent representing the angle between the positive $x$-axis and the line segment from the origin to point $(a,b)$⁠. However, any corresponding selection of $y∗$ that results in $s∗$ in $V˜u$ ensures that $s∗∗$ is not in $V˜u$⁠, since only one point in $V˜u$ is created for each line in the line set. Hence, there is no $y$ that would allow $V˜u≠S$⁠, so $|Au(C,θ,S)|>0$⁠.

APPENDIX B COMPUTATION OF $𝒦^(S_)$

Computation of $𝒦^(S_)$ is simple. Each $𝒦(S_i)$ is the intersection of a set of interior half-planes with one half-plane defined by each edge $S‾i$⁠. For each edge, the plane is divided with the line that intersects the edge. The half-plane that is on the same side of the dividing lines as the interior of the polygon is used. The kernel of $Si$⁠, $𝒦(S_i)$⁠, is the intersection of all these interior half-planes (Shamos, 1975). Intersecting all the individual polygons, $𝒦(S_i)$⁠, produces $𝒦^(S_)$⁠. For any $S_i$ with $n$ edges, $𝒦(S_i)$ can be found as described in $O(nlogn)$ time (Shamos, 1975). An alternative algorithm can find $𝒦(S_i)$ in $O(n)$ time (Lee & Preparata, 1979).

APPENDIX C SIMULATION SHAPE PARAMETERS

This section gives a list of the parameter values used in simulations. All values are rounded to the nearest thousandth. For Shape A, Shape B and Shape C, all the following parameters are shared.

The mean parameters for Shape A, Shape B and Shape C are

Author notes