Skip to content

Instantly share code, notes, and snippets.

@barryrowlingson
Created July 30, 2012 15:39
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save barryrowlingson/3207893 to your computer and use it in GitHub Desktop.
Save barryrowlingson/3207893 to your computer and use it in GitHub Desktop.
Can you fit a rectangle in a polygon?
require(rgeos)
require(sp)
makeRect <- function(origin, l1, l2, angle){
sa = sin(angle)
ca = cos(angle)
x=origin[1] + c(0,l1*ca,l1*ca-l2*sa,-l2*sa,0)
y=origin[2] + c(0,l1*sa,l1*sa+l2*ca,l2*ca,0)
return(SpatialPolygons(list(Polygons(list(Polygon(cbind(x,y))),"Rectangle"))))
return(cbind(origin[1]+x,origin[2]+y))
}
tryRandom <- function(w,h,p,n,all=TRUE){
found=FALSE
pb = bbox(p)
if(all){
foundRects = list()
}
for(i in 1:n){
xr = runif(1,pb[1,1],pb[1,2])
yr = runif(1,pb[2,1],pb[2,2])
angle = runif(1,0,2*pi)
rrect = makeRect(c(xr,yr),w,h,angle)
if(gContains(p,rrect)){
found=TRUE
if(!all){
break
}else{
foundRects[[length(foundRects)+1]]=rrect
}
}
}
if(found){
if(all){
return(foundRects)
}else{
return(rrect)
}
}else{
warning("No rectangle fit")
return(NA)
}
}
# then, using the natural earth lakes data:
# > plot(lakes[125,],asp=1)
# lets see if we can put a 0.1 x 0.003 (degrees!) runway in:
# > rfit = tryRandom(0.1,0.003,lakes[125,],100)
# ignore these warnings, its because we havent set the CRS:
# There were 50 or more warnings (use warnings() to see the first 50)
# if there isn't a warning about not finding any fitting rectangles:
# > invisible(lapply(rfit,plot,add=TRUE,lwd=2,col="#80808080"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment