3

I am trying to calculate the length of shared boundaries between multiple polygons in R. Specifically, I would like to calculate the shared boundary length of all the polygons in Census Place GIS data downloaded from NHGIS from the year 2000. (If you want to obtain the exact data, go to NHGIS.org; Select Data; Under Geographic Levels, select Place; Under years, select 2000; In the table below, select the "GIS Files" tab and download the first file.)

I've looked at this thread, but unfortunately this code throws multiple errors: https://stackoverflow.com/questions/45338384/calculate-the-length-of-shared-boundaries-between-multiple-polygons. I've tried working with the code to modify it but I haven't been able to get it to work. I believe the problem has to do with the complexity of the intersections. One gIntersection often has multiple classes (Points, Lines, Polygons) with different lengths, as seen in the previous link, but the proposed solution does not work.

My question is identical to this unanswered question: Length of the boundary intersection between two polygons in R

    require(rgdal)  
    require(rgeos)
    library(sp)
    library(raster)
    library(spatialEco)
    library(spdep)

    Shapefile <- readOGR(".","US_place_2000")
    Shapefile <- sp.na.omit(Shapefile, margin=1) # NAs are present for some reason?
    Touching_List <- gTouches(Shapefile, byid = TRUE, returnDense=FALSE)
    # also tried Touching_List <- poly2nb(Shapefile, queen = FALSE)

perimeters <- sp::SpatialLinesLengths(as(Shapefile, "SpatialLines"))

    all.length.list <- lapply(1:length(Touching_List), function(from) {
if(is.null(Touching_List[[from]]==FALSE){      
lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE)
      if (class(lines) == "SpatialCollections") {
        list.Lines <- lapply(1:length(Touching_List[[from]]), function(to) {
          line.single <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]][to],])
          if (class(line.single) == "SpatialPoints") {
            # Double the point to create a line
            L1 <- rbind(line.single@coords, line.single@coords)
            rownames(L1) <- letters[1:2]
            Sl1 <- Line(L1)
            Lines.single <- Lines(list(Sl1), ID = as.character(to))
          } else if (class(line.single) == "SpatialLines") {
            Lines.single <- line.single@lines[[1]]
            Lines.single@ID <- as.character(to)
          }
          Lines.single
        })
        lines <- SpatialLines(list.Lines)
      }
      l_lines <- sp::SpatialLinesLengths(lines)
      res <- data.frame(origin = from,
                        perimeter = perimeters[from],
                        touching = Touching_List[[from]],
                        t.length = l_lines,
                        t.pc = 100*l_lines/perimeters[from])
      res
    }
})

    all.length.df <- do.call("rbind", all.length.list)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.