Transformation to sp/sf

Version infromation

In this we expose a list with all the versions of the programms, operative system and the libraries used in this section code, which almost is extracted from sessionInfo().

SO (version): Windows 11 Home (22H2).
R: 4.2.2 (2022-10-31 ucrt).
RStudio: 2022.07.2 (Build 576).
Libraries:
- INLA: INLA_22.05.07
- inlabru: inlabru_2.5.3
- sp: sp_1.5-0
- sf: sf_1.0-8
- ggplot2: ggplot2_3.3.6
- viridis: viridis_0.6.2

There may be cases when you need to work with inla elements as if they were spatial objects (usually sp or sf), this will be particularly appropriate when working with meshes or mesh edges. Whatever the reasons for passing the mesh or mesh edges to spatial objects, there is no built-in function available in library(INLA), so we will proceed with its definition.

Mesh to spatial object

Suppose we start from a set of points given by the XY data frame, on which we apply the function inla.nonconvex.hull to obtain the boundary of our mesh, resulting from the morphological operations of dilation and closure incorporated for the formation of the outer edge.

Get XY data

library(INLA)
library(inlabru)
XY <- readRDS(file="XY.rds")
INLAnonconvex <- inla.nonconvex.hull(XY, convex=1000, resolution=400)
INLAnonconvexMesh <- inla.mesh.2d(boundary=INLAnonconvex, max.edge=5000)

If we plot the mesh ggplot() + inlabru::gg(INLAnonconvexMesh) + theme_void() we obtain the following figure:

Fig. 1: INLA mesh with non convex boundary

Once the mesh is available, to transform the tessellation into a spatial object (sp) we use the information integrated in the mesh object (mesh). Specifically in the rows of the mesh$graph$tv matrix, which store the vertex positions for each i-th triangle. In this way, we simply inscribe the vertex indices for each triangle in the coordinate matrix of the vertices themselves, forming a polygon with each one and finally grouping the result in a SpatialPolygons object, or in a This is synthesized in the following function:

SpatialPolygonDelaunay <- function(mesh){
  return(SpatialPolygons(sapply(1:nrow(mesh$graph$tv), function(i){
    list(Polygons(list(Polygon(mesh$loc[mesh$graph$tv[i,],1:2])),i
    ))
  })
  ))
}

SpatialPolygonDelaunay <- function(mesh){
  return(st_polygon(sapply(1:nrow(mesh$graph$tv), function(i){
    list(mesh$loc[c(mesh$graph$tv[i,], mesh$graph$tv[i,1]),1:2])
  })
  ))
}

Such that the object transformed would be given by:

SpatialPolygonMesh <- SpatialPolygonDelaunay(mesh=INLAnonconvexMesh)

Boundary to spatial object

In this case we intend to transform the edges of the mesh to SpatialPolygons, in such a way that we can account for the gaps present in the mesh. So, the first thing to do is to convert the edges into simply connected polygons (without gaps).

To correctly conform a polygon, the order of succession of the points on the perimeter of the polygon must be indicated, such that the first point is also the last point to close the polygon. This can be easily done because the mesh stores the locations of the vertices that make up the edges as pairs of points in the rows of the mesh$segm$bnd$idx matrix. That is, if the points refer to the same polygon we will have the connection established from the pairs of indexes going through the rows of this matrix, such that it will come to repeat in the first column the position of the initial point to close the polygon.

Therefore, knowing that the sequence of indices for each closed curve implies that the last pair ends at the first point (first pair of indices), we will have that the element of the first column of the next row will be different from the element of the second column of the reference row: mesh$segm$bnd$idx[i+1,1]!=mesh$segm$bnd$idx[i,1]. So, through this condition, we can set the indices of the polygon boundaries that we want to reconstruct as spatial objects (sp or sf).

SpatialPolygonsBoundary <- function(mesh){
  BoundaryIndx <- sapply(1:(nrow(mesh$segm$bnd$idx)-1), function(i){
    return(mesh$segm$bnd$idx[i+1,1]==mesh$segm$bnd$idx[i,2])
  })
  BoundaryIndxLims <- which(!BoundaryIndx)
  BoundaryIndxLims <- c(0, BoundaryIndxLims, nrow(mesh$segm$bnd$idx))
  SpatialPolygonsBoundary <- SpatialPolygons(sapply(2:(length(BoundaryIndxLims)), function(i){
    VertixBoundaryIndex <- matrix(c(mesh$segm$bnd$idx[(BoundaryIndxLims[i-1]+1):BoundaryIndxLims[i],1],
                                    mesh$segm$bnd$idx[BoundaryIndxLims[i],2]), 
                                  ncol=1, byrow=TRUE)
    return(list(Polygons(list(Polygon(mesh$loc[VertixBoundaryIndex,1:2])),i)))
  }))
  return(SpatialPolygonsBoundary)
}

SpatialPolygonsBoundary <- function(mesh){
  BoundaryIndx <- sapply(1:(nrow(mesh$segm$bnd$idx)-1), function(i){
    return(mesh$segm$bnd$idx[i+1,1]==mesh$segm$bnd$idx[i,2])
  })
  BoundaryIndxLims <- which(!BoundaryIndx)
  BoundaryIndxLims <- c(0, BoundaryIndxLims, nrow(mesh$segm$bnd$idx))
  SpatialPolygonsBoundary <- st_polygon(sapply(2:(length(BoundaryIndxLims)), function(i){
    VertixBoundaryIndex <- c(mesh$segm$bnd$idx[(BoundaryIndxLims[i-1]+1),1], 
                             mesh$segm$bnd$idx[(BoundaryIndxLims[i-1]+1):BoundaryIndxLims[i],2])
    return(list(cbind(mesh$loc[VertixBoundaryIndex,1:2])))
  }))
  return(SpatialPolygonsBoundary)
}

Then, once all the polygons have been formed from the edges given by the INLA function, these polygons must be transformed to account for the gaps they may contain and are present. To do this we can either perform a mathematical approach to the properties of the problem, using set theory for ordered sets, or we can approach it operationally to solve it.

Operational reasoning.

To do this according to an operational approach we must consider that we have formed a polygon for each closed curve, for each edge or boundary, so that these curves may contain within them other surfaces given by other curves. In addition, some of these polygons may be associated with the inner edge of the polygons, which would refer to the boundary of the inner hole of the polygon, so that we can have to classify the polygons obtained in two classes or types: full polygons and hollow polygons, depending on whether they refer to the border of a “full” area or a “hollow” area. In addition, it must be taken into account that a filled polygon will contain in the next level of maximum internal intersection, in case of containing in its interior other polygons, a hollow polygon. While for the case of hollow polygons this condition assumes the counterpart; namely, that the next level of maximum internal intersection will be a filled polygon. Such that the level of maximal internal intersectionality with respect to a given polygon we shall understand as the polygon of largest area that is completely contained in the reference polygon. But we can extend this definition according to different sets of mutually independent polygons, where the maximum internal intersectionality for the reference polygon is no longer a single polygon, but the largest polygon of each of these mutually disjoint sets.

On the other hand, if a polygon, bounded by one of these given curves, has a non-zero intersection with others, then this intersection will totally reproduce one of the two polygons.

Subsequently, if analogously to an adjacency or neighborhood matrix, we elaborate an intersectionality matrix, determined according to the intersection operation of a given polygon on the rest of polygons. Such that the matrix will contain the relations between a given polygon with the whole set of polygons according to the intersection operation, operation that will be worth 1 if there is intersection and 0 if there is not. Therefore, for each row or column we will obtain the intersection of a polygon with the rest, where the diagonal of the matrix would manifest the intersection of a polygon with itself, so the main diagonal will be a vector of ones.

From the intersectionality matrix we can conclude that the rows or columns whose sum is unitary (only the element of the main diagonal will be different from 0 for that row or column), will belong to a polygon that will not contain any other polygon, whether it is hollow or full. Therefore, we could extract these polygons as correctly defined polygons that do not require any transformation.

Finally, in case a given row sums to more than one it will imply that it contains at least one hole. So to form the polygons with the holes we can approach the procedure in the following way:

First, for each row we can search for the smallest polygon (area) in the whole set of intersecting polygons, from which a completely ordered set will be constructed according to the intersections of the polygon.
Secondly, the set of polygons can be recomposed as the rearrangement into the smallest number of totally ordered sets. Such that from these totally ordered sets we will know that the largest of the sets refers to a filled polygon, the next to a hollow one, … and so on.
To thirdly, to conform by means of operations of difference and union the different polygons that account for the existing hollows.

Mathematical reasoning.

From a mathematical approach we can use set theory to evaluate the set of operations to be performed to solve the question. To do this, the first thing we must do is to establish the conditions of the problem:

We have a set of $n$ polygons $p_i$, such that the universal set of reference will be $\displaystyle P=\{p_1, …, p_n\}$. These polygons ar esubsets of the two dimensional Euclidian space, $p_i\subset \mathbb{R}^2$.
There are two kinds of polygons: (i) filled polygons $P_{F}=\{p_{fi}\}$ and (ii) hollow polygons $P_H=\{p_{hi}\}$, by jointly reproducing the universal set of polygons $P=P_F\cup P_H$.
There is a function $f$ that stablish the relation between a polygon $p_i$ with its area $a_i$, such that $f:P\longmapsto A$ implies $f(p_i)=a_i$. However, this function is only surjective if we assume that there can be two or more polygons with the same area, while it would be surjective and injective if for each polygon it has a different area.
If the intersection of two polygons is non-zero, then the intersection results in the polygon of smaller area. That is, $p_i\cap p_j=\{\emptyset\veebar p_j\}|f(p_i)>f(p_j)$. In addition, we have as a condition that there is no polygon intersecting with another that has the same area, $(\nexists p_i)(p_i,p_j\in P): f(p_i)=f(p_j)\iff \{p_i\cap p_j\neq \emptyset \wedge i\neq j\}$.
For the set of all non-zero intersections of polygons $(p_i, p_j)$, we will have that the set with the largest area is a filled polygon $f(p_i)>f(p_j):p_i\in P_F$.
Finally, as a weak condition between the non-zero intersection of two filled polygons, we have that there exists, at least, a hollow polygon with greater area than that of the filled polygon of smaller area and that also intersects with them. In other words, if $p_j\subset p_i:p_i,p_j\in P_F\Rightarrow (\exists p_k) p_k\in P_H: p_j\subset p_k\subset p_i$. While the strong condition will be, that if we have two polygons such that one is a subset of another with no filled polygon that can be incorporated as a middle element between the two by the subset operation, then there exists one and only one hollow polygon that lies between the two:

$$(\forall p_i,p_j\in P_F) p_j\subset p_i:\nexists p_k\in P_F, p_j \subset p_k \wedge p_k\subset p_i \iff \exists! p_z\in P_H: p_j\subset p_z \wedge p_z\subset p_i.$$

Given these conditions it is obvious that we can reorder the set $P$ of all the polygons, in such a way that we have the least number of completely ordered sets of higher cardinality whose union reproduces the set { $P$} of reference. That is to say, we have to conform all those completely ordered sets $(P_i,\subset)$, which implies that $(\forall p_i) p_i\in P_i (\exists p_j)p_j\in P_i\rightarrow p_j\subset p_i \vee p_i\subset p_j$. Therefore, we shall construct these completely ordered sets $P_i=\{p_{i1}\subset p_{2i}\subset … \subset p_{ni}\}$, of which we know from conditions (5) and (6) that for any of these sets $P_i$, we will have that the filled polygons will be:

$$ \{ p_{ki},p_{(k-2)i},...,\}\in P_F \Rightarrow p_{(k-2m)i}\in P_F: m\in M\subset\mathbb{N},$$ while the hollow polygons would be: $$ \{ p_{(k-1)i},p_{(k-3)i},...,\}\in P_H \Rightarrow p_{(k-1-2m)i}\in P_H: m'\in M'\subset\mathbb{N}.$$

The code to perform the pertinent transformations and obtain the polygons with the corresponding holes is as follows:

CompleteOrderedSetsFromSpatialPolygonsClass <- function(SPDINLAnonconvexBoundary){
  SPDINLABoundaryList <- list(Polygons=list(), 
                              DFIDArea=data.frame(matrix(ncol=2,nrow=0, dimnames=list(NULL, c("ID", "Area"))))
  )
  for(i in 1:length(SPDINLAnonconvexBoundary)){
    SPDINLABoundaryList$Polygons[[i]] <- SPDINLAnonconvexBoundary[i]
    SPDINLABoundaryList$DFIDArea[i,] <- c(i, SPDINLAnonconvexBoundary[i]@polygons[[1]]@area)
  }
  
  IndexIDV <- SPDINLABoundaryList$DFIDArea$ID
  
  ConnectivityMatrix <- matrix(ncol=length(SPDINLAnonconvexBoundary), nrow=length(SPDINLAnonconvexBoundary))
  for(i in 1:length(SPDINLAnonconvexBoundary)){
    for(j in i:length(SPDINLAnonconvexBoundary)){
      ConnectivityMatrix[i,j] <- ConnectivityMatrix[j,i] <- as.numeric(gIntersects(SPDINLAnonconvexBoundary[i],
                                                                                   SPDINLAnonconvexBoundary[j]))
    }
  }
  
  IndexIsolatedPolygonFill <- which(apply(ConnectivityMatrix, FUN=sum, MARGIN=1)==1) 

  IndexIDV <- if(length(IndexIsolatedPolygonFill)==0){IndexIDV} else{IndexIDV[-which(apply(ConnectivityMatrix, FUN=sum, MARGIN=1)==1)]}
  
  CompleteOrderedSets <- list(CompleteOrderedSets=list(), DFCompleteOrderedSets=list(), 
                              CompleteOrderedIsolatedSets=list(), DFCompleteOrderedIsolatedSets=list())
  
  if(length(IndexIsolatedPolygonFill)!=0){
    for(i in 1:length(IndexIsolatedPolygonFill)){
      CompleteOrderedSets$CompleteOrderedIsolatedSets[[i]] <- SPDINLABoundaryList$Polygons[[IndexIsolatedPolygonFill[i]]]
      CompleteOrderedSets$DFCompleteOrderedIsolatedSets[[i]] <- 
        data.frame(ID=SPDINLABoundaryList$Polygons[[IndexIsolatedPolygonFill[i]]]@polygons[[1]]@ID,
                   Area=SPDINLABoundaryList$Polygons[[IndexIsolatedPolygonFill[i]]]@polygons[[1]]@area,
                   HOLE=FALSE)
    }
  }
  
  if(length(IndexIsolatedPolygonFill)!=length(IndexIDV)){
    n_step <- 1
    while(length(IndexIDV)!=0){
      IDminArea <- SPDINLABoundaryList$DFIDArea[IndexIDV,"ID"][which.min(SPDINLABoundaryList$DFIDArea[IndexIDV,"Area"])] 
      IndxCOS <- which(as.logical(ConnectivityMatrix[IDminArea,]))
      DFIDAreaV <- SPDINLABoundaryList$DFIDArea[IndxCOS,]
      CompleteOrderedSets$CompleteOrderedSets[[n_step]] <- 
        SpatialPolygons(sapply(DFIDAreaV[order(DFIDAreaV$Area, decreasing=TRUE),"ID"],
                               function(i){
                                 list(SPDINLABoundaryList$Polygons[[i]]@polygons[[1]])
                               })
        )
      CompleteOrderedSets$DFCompleteOrderedSets[[n_step]] <- 
        data.frame(ID=DFIDAreaV[order(DFIDAreaV$Area, decreasing=TRUE),"ID"],
                   Area=DFIDAreaV[order(DFIDAreaV$Area, decreasing=TRUE),"Area"],
                   HOLE=floor(1:length(IndxCOS)/2)==1:length(IndxCOS)/2
        )
      
      IndexIDV <- setdiff(IndexIDV,IndxCOS)
      n_step <- n_step+1
    } 
    SpatialPolygonsNames <- c()
    for(i in 1:length(CompleteOrderedSets$CompleteOrderedSets)){
      if(length(CompleteOrderedSets$CompleteOrderedSets[[i]])/2!=floor(length(CompleteOrderedSets$CompleteOrderedSets[[i]])/2)){
        SpatialPolygonsNames <- 
          c(SpatialPolygonsNames, 
            paste0("SPDINLA", 
                   CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[length(CompleteOrderedSets$CompleteOrderedSets[[i]])]))
        
        assign(x=paste0("SPDINLA",
                        CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[length(CompleteOrderedSets$CompleteOrderedSets[[i]])]),
               value=SpatialPolygons(list(
                 CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[length(CompleteOrderedSets$CompleteOrderedSets[[i]])]]
               )))
      }
      
      for(j in 1:floor(length(CompleteOrderedSets$CompleteOrderedSets[[i]])/2)){
        indx <- 1:2+2*(j-1)
        if(exists(paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]))){
          assign(x=paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]),
                 value=gDifference(eval(parse(text=paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]))), 
                                   SpatialPolygons(list(CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[indx[2]]])),
                                   id=paste0(CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]])
                 )
          )
        } else{
          SpatialPolygonsNames <- c(SpatialPolygonsNames, paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]))
          assign(x=paste0("SPDINLA", CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]]),
                 value=gDifference(SpatialPolygons(list(CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[indx[1]]])), 
                                   SpatialPolygons(list(CompleteOrderedSets$CompleteOrderedSets[[i]]@polygons[[indx[2]]])),
                                   id=paste0(CompleteOrderedSets$DFCompleteOrderedSets[[i]]$ID[indx[1]])
                 )
          )
        }}
      
      if(length(IndexIsolatedPolygonFill)!=0){
        SpatialPolygonRefined <- SpatialPolygons(c(lapply(1:length(SpatialPolygonsNames), function(i){
          return(eval(parse(text=paste0(SpatialPolygonsNames[i], "@polygons[[1]]"))))
        }), lapply(1:length(CompleteOrderedSets$CompleteOrderedIsolatedSets), function(i){
          return(CompleteOrderedSets$CompleteOrderedIsolatedSets[[i]]@polygons[[1]])
        })))        
      } else{
        SpatialPolygonRefined <- SpatialPolygons(lapply(1:length(SpatialPolygonsNames), function(i){
          return(eval(parse(text=paste0(SpatialPolygonsNames[i], "@polygons[[1]]"))))
        }))
      }

    }
  } else{
    SpatialPolygonRefined <- SpatialPolygons(lapply(1:length(CompleteOrderedSets$CompleteOrderedIsolatedSets), function(i){
      return(CompleteOrderedSets$CompleteOrderedIsolatedSets[[i]]@polygons[[1]])
    }))
  }
  
  SpatialPolygonRefined@proj4string <- CRS("+proj=utm")
  CompleteOrderedSetsAndPolygons <- list(CompleteOrderedSets=CompleteOrderedSets,
                                         SpatialPolygonRefined=SpatialPolygonRefined)
  
  return(CompleteOrderedSetsAndPolygons)
}

Thanks to this function we have all the polygons according to the closed curves that determine them. This allows us to perform operations on them, so that we are taking into account the existing gaps.

Therefore, following the example, we would have to apply the first function SpatialPolygonsBoundary on the result of the mesh INLAnonconvexMesh, and then apply the function just defined.

SPDINLAnonconvexBoundary <- SpatialPolygonsBoundary(mesh = INLAnonconvexMesh)
PolygonsRefined <- CompleteOrderedSetsFromSpatialPolygonsClass(SPDINLAnonconvexBoundary = SPDINLAnonconvexBoundary)

In order to visualize that we now have all the polygons available independent of each other, we can plot them in different colors:

col <- viridis::turbo(length(PolygonsRefined$SpatialPolygonRefined))
ggplotPolygons <- ggplot()
for(i in 1:length(PolygonsRefined$SpatialPolygonRefined)){
    ggplotPolygons <- ggplotPolygons + geom_sf(data=st_as_sf(PolygonsRefined$SpatialPolygonRefined[i]), fill=col[i], alpha=0.5)
}

So, if we represent the polygons as ggplotPolygons + theme_void() we will obtain the following figure:

Fig. 2: Refined polygons

Example with nested polygons.