[R-sig-Geo] RE : RE : determining suitable lakes for landing

[Bastien.Ferland-Raymond at mrnf.gouv.qc.ca]

So, I've applied Barry's approach to one of my test territory. Barry nicely qualified is approach as brute-force, so I've made my computer sweat a little.  I've tried 20000 iterations per polygon.  It was surprisingly not that hard.  First, ~90% of my lakes were remove through simple Pythagorean analyses, leaving 162 lakes.  For those lakes, some where big and went through Barry's algorithm very fast. Finally, I speed up the process using parallel processing to take advantage of my 8 cores processor.  All that together, we are talking a 20 minutes analysis, which is really viable in my context.  Here is my code :

#############
library(rgdal)
lacs <- readOGR("M:\\carte R\\en mtm","lacs_amerrissables")

##### begin fonction
ameri <- function(nom, lacs.shape,iter) {
  require(rgeos); require(sp)

  trop.petit <- function(lac) 1500>sqrt(sum((bbox(lac)[,2]-bbox(lac)[,1])^2))

  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")), proj4string=CRS(proj4string(lac))))
    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){
#    dans.pol=F
#    while(dans.pol==F){
      xr = runif(1,pb[1,1],pb[1,2])
      yr = runif(1,pb[2,1],pb[2,2])
#     dans.pol <- gIntersects(SpatialPoints(cbind(xr,yr), proj4string=CRS(proj4string(lac))),lac)
#    }
      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(NULL)
    }
  }

 lac <- lacs.shape[lacs.shape at data$NOUNIQUE==nom,]
 if(trop.petit(lac)) {
   land <- "non"
  } else {
   nul <- !is.null(tryRandom(50,1500,lac,iter, all=F))
   land <- ifelse(nul, "oui", "p-e")
  }
 land
}
#### end function

ameri(lacs at data$NOUNIQUE[1666], lacs, iter=30000)  # for one lake

library(parallel)
cl=makeCluster(7)
system.time(is.ameri <- parSapply(cl,lacs at data$NOUNIQUE,ameri, lacs.shape=lacs, iter=20000))    ## for multiple lakes
stopCluster(cl)

###############

The code gave pretty good results, however not perfect.  27 lakes gave different results from the initial Arcinfo code.  Following manually looking at those lakes, I've noticed that 40% were missed by Barry's algorithm (false negative).  However, that also meant that the ArcInfo was doing 60% of false-positive, which makes Barry's code better :)  Increasing the number of simulations to 30000 helped finding some of those.

One thing I'm thinking right now that could help the process is to make sure every rectangle origin falls within the lake instead of within the extent of the lake.  This may help make sure it tests only possible rectangle, however the speed may not be that different as it can be kind of hard to find points in the lake in the case of very long narrow lake.

Finally, I'll still have to work on the island problem.  I'll look into rgeos hole managing capability.

Concerning Karl Ove Hufthammer idea.  It's seems nice, however I'm afraid I won't be able to put the time to understand those algorithms and implement them in R.

Thanks again for your help,

Bastien