Testing if geometries have at least one boundary point in common regardless of interior points (not like `gTouches`)

Question

the gTouches function in the rgeos package tests whether "geometries have at least one boundary point in common, but no interior points". I am looking for a way to test whether "geometries have at least one boundary point in common" without the criteria related to interior points.

Here is the basic setup: I have two shapefiles that are mostly embedded in each other. I want to find the polygons in the file with the smaller areas that are at the border of the larger areas. Here is a graph to describe what I am trying to do:

plot(map2, col=NA, border='black', lwd=0.4)
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)

The figure shows census blocks in Staten Island, NY. The green highlighting in one of the larger areas illustrates the blocks I want to identify. Only those that share or cross a border of the larger areas (thick lines). Not the blocks that are in the middle of the larger areas. I tried to do this with with gTouches(map2,map1, byid=TRUE) and other function in the rgeos package but without success. gTouches only returns FALSE probably because the criteria is that "geometries have at least one boundary point in common, but no interior points". Basically, I am looking for a function that tests whether "geometries have at least one boundary point in common" regardless of the interior.

A follow-up question is whether I can get the length of the mutual border?

Data: You can download the two maps here and here. Both are rds files so you can open them like this:

library('rgdal')
library('rgeos')
library('sp')
map1 = readRDS('map1.rds')
map2 = readRDS('map2.rds')

Answer 1

You can use a combo of gIntersects() (to find all little polygons that intersect any part of the school district) and gContainsProperly() (to find all little polygons that are fully contained within and not intersecting the boundary of the school district). Then simply combine the two resulting logical matrices to identify the polygons you're after.

## Identify polygons that intersect but aren't fully contained within the
## school district whose polygon is given by SD = map2[13,]
SD <- map2[13,]
ii <- gIntersects(SD, map1, byid=TRUE) &
      !gContainsProperly(SD, map1, byid=TRUE)
ii <- apply(ii, 1, any)  ## Handy construct if both layers contain >1 polygon

## Plot that area, to show that results are  correct
plot(SD, col=NA, border='black')                  ## Establish plotted area
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)
plot(map1[ii,], col="lightgreen", add=TRUE)
plot(SD, col=NA, border='black', lwd=2, add=TRUE) ## Put SD boundary on top 

EDIT :

That's not quite right, however. As can be seen in the map above, many small polygons along the SW and SE interiors of the school district which should have been identified have not been. Outcomes like this occur pretty frequently with rgeos operations, and arise from tiny misregistrations of the pair of layers (or of their intermediate representations by the GEOS engine).

The solution is to use gBuffer() to buffer out one of the layers by a small amount before performing the topological queries. Here, the coordinates are in meters, and a bit of trial and error showed that a 20-meter buffer turns out to be mostly sufficient to fix the problem:

## Expand every polygon in map1 by a 20-meter wide buffer 
map1_buff <- gBuffer(map1, byid=TRUE, width=20)

## and then use the buffered version of map1 in the topological queries
ii <- gIntersects(SD, map1_buff, byid=TRUE) &
      !gContainsProperly(SD, map1_buff, byid=TRUE)
ii <- apply(ii, 1, any)  ## Handy construct if both layers contain >1 polygon

## Plot that area, to show that results are  correct
plot(SD, col=NA, border='black')                  ## Establish plotted area
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)
plot(map1[ii,], col="lightgreen", add=TRUE)
plot(SD, col=NA, border='black', lwd=2, add=TRUE)

This still misses a couple of polygons along the coast, but at some point a complete solution may have to involve getting a pair of maps that match up better in their level of detail. If the buffer size is made much larger, this analysis will start to produce false positives, picking up, for example, a few of the truly interior polygons in the NW corner of the school district.