Various Functions

Here we will present different functions that have been developed for implementation in automated processes or that may be of interest.

  • Transformation to sp/sf

    Here we define some functions related to the transformations of inla.mesh objects to sp and sf.

Subsections of Various Functions

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
INLA_mesh_with_non_convex_boundary.png

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.
Mathematical reasoning.

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
Polygons_with_holes
Example with nested polygons.