sl.polygon.mergehole <-
function(poly,poly.hole,connect.maxstep=Inf,checkcross.input=TRUE,checkcross.output=TRUE,split.poly=FALSE,checkfail.maskfrac=0.01,split.costexp=2,split.costborderbygcdist=FALSE) {
if (any(!(c("lon","lat") %in% names(poly.hole)))) {
stop("handling more than one hole at a time has not yet been implemented")
}
o.N = length(poly$lon)
dupl = which(poly$lon == c(poly$lon[o.N],poly$lon[1:(o.N-1)])
& poly$lat == c(poly$lat[o.N],poly$lat[1:(o.N-1)]))
if (length(dupl) > 0) {
warning(paste0("removing ",length(dupl)," duplicate points from 'poly'"))
poly$lon = poly$lon[-dupl]
poly$lat = poly$lat[-dupl]
o.N = length(poly$lon)
}
if (checkcross.input) {
intersect.res = sl.intersect(poly$lon,poly$lat,close1=TRUE,return.on.firsthit=TRUE)
if (intersect.res$anylines.intersect) {browser();stop("'poly' contains intersecting segments")}
}
o.lon.ext = c(poly$lon[o.N], poly$lon, poly$lon[1])
o.lat.ext = c(poly$lat[o.N], poly$lat, poly$lat[1])
o.angles = numeric(o.N)
for (i in 1:o.N) {
o.angles[i] = sl.angle(o.lon.ext[i:(i+2)],o.lat.ext[i:(i+2)],left=TRUE)
}
if (sum(o.angles) > o.N*pi) {
poly$lon = poly$lon[o.N:1]
poly$lat = poly$lat[o.N:1]
o.lon.ext = c(poly$lon[o.N], poly$lon, poly$lon[1])
o.lat.ext = c(poly$lat[o.N], poly$lat, poly$lat[1])
o.angles = 2*pi - o.angles[o.N:1]
}
i.N = length(poly.hole$lon)
dupl = which(poly.hole$lon == c(poly.hole$lon[i.N],poly.hole$lon[1:(i.N-1)])
& poly.hole$lat == c(poly.hole$lat[i.N],poly.hole$lat[1:(i.N-1)]))
if (length(dupl) > 0) {
warning(paste0("removing ",length(dupl)," duplicate points from 'poly.hole'"))
poly.hole$lon = poly.hole$lon[-dupl]
poly.hole$lat = poly.hole$lat[-dupl]
i.N = length(poly.hole$lon)
}
if (checkcross.input) {
intersect.res = sl.intersect(poly.hole$lon,poly.hole$lat,close1=TRUE,return.on.firsthit=TRUE)
if (intersect.res$anylines.intersect) {stop("'poly.hole' contains intersecting segments")}
}
i.lon.ext = c(poly.hole$lon[i.N], poly.hole$lon, poly.hole$lon[1])
i.lat.ext = c(poly.hole$lat[i.N], poly.hole$lat, poly.hole$lat[1])
i.angles = numeric(i.N)
for (i in 1:i.N) {
i.angles[i] = sl.angle(i.lon.ext[i:(i+2)],i.lat.ext[i:(i+2)],left=TRUE)
}
if (sum(i.angles) < i.N*pi) {
poly.hole$lon = poly.hole$lon[i.N:1]
poly.hole$lat = poly.hole$lat[i.N:1]
i.lon.ext = c(poly.hole$lon[i.N], poly.hole$lon, poly.hole$lon[1])
i.lat.ext = c(poly.hole$lat[i.N], poly.hole$lat, poly.hole$lat[1])
i.angles = 2*pi - i.angles[i.N:1]
}
dist.res = sl.finddist(poly$lon,poly$lat,poly.hole$lon,poly.hole$lat,return.alldist=TRUE)
o.ind = dist.res$ind[1]
i.ind = dist.res$ind[2]
connect.dist = dist.res$dist
o.dist = NULL
checkfail.maskN.o = max(floor((floor(checkfail.maskfrac * o.N) - 1) / 2), 0)
checkfail.maskN.i = max(floor((floor(checkfail.maskfrac * i.N) - 1) / 2), 0)
repeat {
o.inds = NULL
if (o.ind < o.N) {o.inds = (o.ind+1):o.N}
if (o.ind > 1) {o.inds = c(o.inds,1:(o.ind-1))}
i.inds = NULL
if (i.ind < i.N) {i.inds = (i.ind+1):i.N}
if (i.ind > 1) {i.inds = c(i.inds,1:(i.ind-1))}
if (connect.dist == 0) {break}
i.ang = sl.angle(c(i.lon.ext[c(i.ind,i.ind+1)],poly$lon[o.ind]),
c(i.lat.ext[c(i.ind,i.ind+1)],poly$lat[o.ind]),left=TRUE)
o.ang = sl.angle(c(o.lon.ext[c(o.ind,o.ind+1)],poly.hole$lon[i.ind]),
c(o.lat.ext[c(o.ind,o.ind+1)],poly.hole$lat[i.ind]),left=TRUE)
angles.OK = (i.ang >= 0 && i.ang <= i.angles[i.ind] && o.ang >= 0 && o.ang <= o.angles[o.ind])
nocross = TRUE
if (angles.OK && checkcross.output) {
if (is.null(o.dist)) {
o.dist = sl.gc.dist(o.lon.ext[2:(o.N+2)], o.lat.ext[2:(o.N+2)], sequential=TRUE)
i.dist = sl.gc.dist(i.lon.ext[2:(i.N+2)], i.lat.ext[2:(i.N+2)], sequential=TRUE)
}
connect.lon = c(poly.hole$lon[i.ind],poly$lon[o.ind])
connect.lat = c(poly.hole$lat[i.ind],poly$lat[o.ind])
intersect.res = sl.intersect(connect.lon, connect.lat, poly.hole$lon[i.inds], poly.hole$lat[i.inds],
dist1 = connect.dist, dist2 = i.dist[i.inds[-length(i.inds)]],
return.on.firsthit = TRUE)
if (!intersect.res$anylines.intersect) {
intersect.res = sl.intersect(connect.lon, connect.lat, poly$lon[o.inds], poly$lat[o.inds],
dist1 = connect.dist, dist2 = o.dist[o.inds[-length(o.inds)]],
return.on.firsthit = TRUE)
}
nocross = !intersect.res$anylines.intersect
}
if (angles.OK && nocross) {break}
if (checkfail.maskN.o > 0) {
o.ind = (((o.ind-checkfail.maskN.o-1):(o.ind+checkfail.maskN.o-1)) %% o.N) + 1
}
if (checkfail.maskN.i > 0) {
i.ind = (((i.ind-checkfail.maskN.i-1):(i.ind+checkfail.maskN.i-1)) %% i.N) + 1
}
dist.res$dist.all[o.ind,i.ind] = NA
connect.dist = min(dist.res$dist.all,na.rm=TRUE)
oi.inds = which(dist.res$dist.all == connect.dist, arr.ind = TRUE)[1,]
o.ind = oi.inds[1]
i.ind = oi.inds[2]
}
if (is.function(connect.maxstep)) {
if (is.null(o.dist)) {
o.dist = sl.gc.dist(o.lon.ext[2:(o.N+2)], o.lat.ext[2:(o.N+2)], sequential=TRUE)
i.dist = sl.gc.dist(i.lon.ext[2:(i.N+2)], i.lat.ext[2:(i.N+2)], sequential=TRUE)
}
connect.maxstep = connect.maxstep(c(o.dist,i.dist))
}
if (connect.dist > connect.maxstep) {
np = ceiling(connect.dist/connect.maxstep) - 1
connect.insert = sl.fillequidist(c(poly$lon[o.ind],poly.hole$lon[i.ind]),
c(poly$lat[o.ind],poly.hole$lat[i.ind]),np=np,method="gc")
connect.insert.inv = list(lon=connect.insert$lon[np:1],lat=connect.insert$lat[np:1])
} else {
connect.insert = list(lon=NULL,lat=NULL)
connect.insert.inv = list(lon=NULL,lat=NULL)
}
if (!split.poly) {
lon = c(poly$lon[c(o.ind,o.inds,o.ind)], connect.insert$lon,
poly.hole$lon[c(i.ind,i.inds,i.ind)], connect.insert.inv$lon)
lat = c(poly$lat[c(o.ind,o.inds,o.ind)], connect.insert$lat,
poly.hole$lat[c(i.ind,i.inds,i.ind)], connect.insert.inv$lat)
return(list(lon=lon, lat=lat))
}
o.ind.prev = o.ind
i.ind.prev = i.ind
connect.dist.prev = connect.dist
if (checkcross.output) {
connect.lon.prev = connect.lon
connect.lat.prev = connect.lat
}
connect.insert.prev = connect.insert
connect.insert.inv.prev = connect.insert.inv
if (split.costborderbygcdist) {
o.dist.cumsum = cumsum(o.dist)
i.dist.cumsum = cumsum(i.dist)
} else {
o.dist.cumsum = 1:o.N
i.dist.cumsum = 1:o.N
}
o.dist.total = tail(o.dist.cumsum,1)
o.dist.cumsum = (o.dist.cumsum - o.dist.cumsum[o.ind.prev])[-o.ind.prev]
o.dist.min = pmin(o.dist.total-o.dist.cumsum, abs(o.dist.cumsum), o.dist.total+o.dist.cumsum)
i.dist.total = tail(i.dist.cumsum,1)
i.dist.cumsum = (i.dist.cumsum - i.dist.cumsum[i.ind.prev])[-i.ind.prev]
i.dist.min = pmin(i.dist.total-i.dist.cumsum, abs(i.dist.cumsum), i.dist.total+i.dist.cumsum)
split.costfun = dist.res$dist.all[-o.ind.prev,-i.ind.prev] / (pmin(matrix(o.dist.min,nrow=o.N-1,ncol=i.N-1),
matrix(i.dist.min,nrow=o.N-1,ncol=i.N-1,byrow=TRUE)))^split.costexp
oi.inds = which(split.costfun == min(split.costfun,na.rm=TRUE), arr.ind = TRUE)[1,]
o.ind = oi.inds[1]
if (o.ind >= o.ind.prev) {o.ind = o.ind + 1}
i.ind = oi.inds[2]
if (i.ind >= i.ind.prev) {i.ind = i.ind + 1}
connect.dist = dist.res$dist.all[o.ind,i.ind]
repeat {
o.inds = NULL
if (o.ind < o.N) {o.inds = (o.ind+1):o.N}
if (o.ind > 1) {o.inds = c(o.inds,1:(o.ind-1))}
i.inds = NULL
if (i.ind < i.N) {i.inds = (i.ind+1):i.N}
if (i.ind > 1) {i.inds = c(i.inds,1:(i.ind-1))}
if (connect.dist == 0) {break}
i.ang = sl.angle(c(i.lon.ext[c(i.ind,i.ind+1)],poly$lon[o.ind]),
c(i.lat.ext[c(i.ind,i.ind+1)],poly$lat[o.ind]),left=TRUE)
o.ang = sl.angle(c(o.lon.ext[c(o.ind,o.ind+1)],poly.hole$lon[i.ind]),
c(o.lat.ext[c(o.ind,o.ind+1)],poly.hole$lat[i.ind]),left=TRUE)
angles.OK = (i.ang >= 0 && i.ang <= i.angles[i.ind] && o.ang >= 0 && o.ang <= o.angles[o.ind])
nocross = TRUE
if (angles.OK && checkcross.output) {
connect.lon = c(poly.hole$lon[i.ind],poly$lon[o.ind])
connect.lat = c(poly.hole$lat[i.ind],poly$lat[o.ind])
intersect.res = sl.intersect(connect.lon, connect.lat, poly.hole$lon[i.inds], poly.hole$lat[i.inds],
dist1 = connect.dist, dist2 = i.dist[i.inds[-length(i.inds)]],
return.on.firsthit = TRUE)
if (!intersect.res$anylines.intersect) {
intersect.res = sl.intersect(connect.lon, connect.lat, poly$lon[o.inds], poly$lat[o.inds],
dist1 = connect.dist, dist2 = o.dist[o.inds[-length(o.inds)]],
return.on.firsthit = TRUE)
}
if (!intersect.res$anylines.intersect) {
intersect.res = sl.intersect(connect.lon, connect.lat, connect.lon.prev, connect.lat.prev,
dist1 = connect.dist, dist2 = connect.dist.prev,
return.on.firsthit = TRUE)
}
nocross = !intersect.res$anylines.intersect
}
if (angles.OK && nocross) {break}
o.ind = oi.inds[1]
i.ind = oi.inds[2]
if (checkfail.maskN.o > 0) {
o.ind = (((o.ind-checkfail.maskN.o-1):(o.ind+checkfail.maskN.o-1)) %% o.N) + 1
}
if (checkfail.maskN.i > 0) {
i.ind = (((i.ind-checkfail.maskN.i-1):(i.ind+checkfail.maskN.i-1)) %% i.N) + 1
}
split.costfun[o.ind,i.ind] = NA
oi.inds = which(split.costfun == min(split.costfun,na.rm=TRUE), arr.ind = TRUE)[1,]
o.ind = oi.inds[1]
if (o.ind >= o.ind.prev) {o.ind = o.ind + 1}
i.ind = oi.inds[2]
if (i.ind >= i.ind.prev) {i.ind = i.ind + 1}
connect.dist = dist.res$dist.all[o.ind,i.ind]
}
if (connect.dist > connect.maxstep) {
np = ceiling(connect.dist/connect.maxstep) - 1
connect.insert = sl.fillequidist(c(poly$lon[o.ind],poly.hole$lon[i.ind]),
c(poly$lat[o.ind],poly.hole$lat[i.ind]),np=np,method="gc")
connect.insert.inv = list(lon=connect.insert$lon[np:1],lat=connect.insert$lat[np:1])
} else {
connect.insert = list(lon=NULL,lat=NULL)
connect.insert.inv = list(lon=NULL,lat=NULL)
}
o.inds.cut = which(o.inds == o.ind.prev)
o.inds.1 = NULL
if (o.inds.cut < (o.N-1)) {o.inds.1 = o.inds[(o.inds.cut+1):(o.N-1)]}
o.inds.2 = NULL
if (o.inds.cut > 1) {o.inds.2 = o.inds[1:(o.inds.cut-1)]}
i.inds.cut = which(i.inds == i.ind.prev)
i.inds.1 = NULL
if (i.inds.cut > 1) {i.inds.1 = i.inds[1:(i.inds.cut-1)]}
i.inds.2 = NULL
if (i.inds.cut < (i.N-1)) {i.inds.2 = i.inds[(i.inds.cut+1):(i.N-1)]}
lon1 = c(poly$lon[c(o.ind.prev,o.inds.1,o.ind)], connect.insert$lon,
poly.hole$lon[c(i.ind,i.inds.1,i.ind.prev)], connect.insert.inv.prev$lon)
lat1 = c(poly$lat[c(o.ind.prev,o.inds.1,o.ind)], connect.insert$lat,
poly.hole$lat[c(i.ind,i.inds.1,i.ind.prev)], connect.insert.inv.prev$lat)
lon2 = c(poly$lon[c(o.ind,o.inds.2,o.ind.prev)], connect.insert.prev$lon,
poly.hole$lon[c(i.ind.prev,i.inds.2,i.ind)], connect.insert.inv$lon)
lat2 = c(poly$lat[c(o.ind,o.inds.2,o.ind.prev)], connect.insert.prev$lat,
poly.hole$lat[c(i.ind.prev,i.inds.2,i.ind)], connect.insert.inv$lat)
return(list(polygon1=list(lon=lon1, lat=lat1), polygon2=list(lon=lon2, lat=lat2)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.