Dear R gurus,
I'm facing a challenge that I have a good feeling you could help me with. So far I'm stuck right at the beginning so I don't have any example to show you. Here is what I'm trying to do: I have a polygon shape of lakes for which I have to be able to automatically determine if they are suitable (large enough) for a plane to land. To be suitable, we must be able to fit a rectangle of dimension 50m x 1500m in the lake without it touching any lake contour. Of course, the rectangle can be in any orientation in the lake, which makes the test a little harder. I can easily eliminate all the small lakes by using the extent and Pythagorean. For exemple: 1500 > sqrt(sum((bbox(lake)[,2] - bbox(lake)[,1])^2)) # measure the diagonal of the bbox of the lake (mtm projection) All those lake are definitely too small for a plane to land, however not all the lakes that are bigger than that will be large enough. Anybody have an idea how to place the rectangle and discriminate those lakes? Thanks for your help! Bastien Bastien Ferland-Raymond, M.Sc. Stat., M.Sc. Biol. Division des orientations et projets spéciaux Direction des inventaires forestiers Ministère des Ressources naturelles et de la Faune _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
On Mon, Jul 30, 2012 at 1:11 PM,
[hidden email] <[hidden email]> wrote: > Dear R gurus, > > I'm facing a challenge that I have a good feeling you could help me with. So far I'm stuck right at the beginning so I don't have any example to show you. Here is what I'm trying to do: > > I have a polygon shape of lakes for which I have to be able to automatically determine if they are suitable (large enough) for a plane to land. To be suitable, we must be able to fit a rectangle of dimension 50m x 1500m in the lake without it touching any lake contour. Of course, the rectangle can be in any orientation in the lake, which makes the test a little harder. > > I can easily eliminate all the small lakes by using the extent and Pythagorean. For exemple: > > 1500 > sqrt(sum((bbox(lake)[,2] - bbox(lake)[,1])^2)) # measure the diagonal of the bbox of the lake (mtm projection) > > All those lake are definitely too small for a plane to land, however not all the lakes that are bigger than that will be large enough. > > Anybody have an idea how to place the rectangle and discriminate those lakes? > > Thanks for your help! Fun problem... It seems there are algorithms for finding the maximum area rectangle [http://www.sciencedirect.com/science/article/pii/S0096300312003207] but in your case the max area rectangle could be 49m x 10000m, and the lake could have a 50mx1500m finger poking out of it. The max area rectangle would then not be good enough to land in, but the lake would actually be okay. But if the max area rectangle is > 50m x > 1500m then you've got one you can land in, but its not a necessary condition.... Some kind of heuristic approach of dropping 50mx1500m rectangles on the lake area and testing might work - you'd have to randomise the location and orientation and vary the randomness by how much the rectangle overlaps, some kind of simulated annealing approach perhaps. Again, its not guaranteed to work due to its heuristic approach - eventually your simulated annealing has to stop and it might not have found a solution... In the pathological case I mentioned above it might have trouble... I think that if a polygon contains an NxM rectangle then that rectangle can be shifted to touch two points of the polygon, so maybe you can reduce the problem to rectangles that touch at two points - that might help any heuristic algorithms... I wonder if everyone on R-sig-geo is now sliding bus tickets around on pieces of paper like me... Barry _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
In reply to this post by Bastien.Ferland-Raymond
On Mon, Jul 30, 2012 at 1:52 PM, Rowlingson, Barry
<[hidden email]> wrote: >> Anybody have an idea how to place the rectangle and discriminate those lakes? To save pasting masses of code in here, I've put a random rectangle tester in a gist: https://gist.github.com/3207893 that does a stupid brute-force approach of trying random NxM rectangles at random angles all over the polygon until it finds one or gives up. I reckon any polygon that doesn't get a solution after 5000 trials is probably so horrendous that no pilot would want to land there anyway... Barry _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Thanks Barry for your interest in my problem. It is indeed an interesting problem at the same time as a complex one.
I've tried your code and it works for two lakes so far. It's seems like a very good base for what I want to do. I'll modify it a little to speed it up. First by removing the very small lakes, than by changing the loop "for" by a "while" as I just need to find if I can fit or not a rectangle, I don't need the rectangle(s). Here is a example of shape I'm going to apply it to: https://www.dropbox.com/s/24eg0eetm0o6dth/lacs_amerrissables.zip There is already a field telling if you could land of not ("AMERISS") that was found by some old ARCgis coding (ARCinfo and aml, not sure exactly how). It was a consultant that did this code for us a couple of years ago, now we want to transfer to R, explaining why I have this problem now. The code was also using iteration, however not random, on rasterized lakes. I think it would be interesting to compare the result from the 2 techniques. I'm going to try that this afternoon, and let you know the result. Thanks again! Bastien -----Message d'origine----- De : [hidden email] [mailto:[hidden email]] De la part de Barry Rowlingson Envoyé : 30 juillet 2012 11:41 À : Ferland-Raymond, Bastien (DIF) Cc : [hidden email] Objet : Re: [R-sig-Geo] determining suitable lakes for landing On Mon, Jul 30, 2012 at 1:52 PM, Rowlingson, Barry <[hidden email]> wrote: >> Anybody have an idea how to place the rectangle and discriminate those lakes? To save pasting masses of code in here, I've put a random rectangle tester in a gist: https://gist.github.com/3207893 that does a stupid brute-force approach of trying random NxM rectangles at random angles all over the polygon until it finds one or gives up. I reckon any polygon that doesn't get a solution after 5000 trials is probably so horrendous that no pilot would want to land there anyway... Barry _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
In reply to this post by Barry Rowlingson
On Mon, Jul 30, 2012 at 5:16 PM,
[hidden email] <[hidden email]> wrote: > Thanks Barry for your interest in my problem. It is indeed an interesting problem at the same time as a complex one. > > I've tried your code and it works for two lakes so far. It's seems like a very good base for what I want to do. I'll modify it a little to speed it up. First by removing the very small lakes, than by changing the loop "for" by a "while" as I just need to find if I can fit or not a rectangle, I don't need the rectangle(s). There's an option to my function that you can use to tell it to stop if it finds it can fit a rectangle - if you pass all=FALSE it returns as soon as it finds one, otherwise it returns NA. I'm not sure how my method will work if there are islands (I assume you dont want islands in the middle of your landing strip!). It uses rgeos:gContains for the test, so I think as long as the island is properly coded as an island it wont consider rectangles with islands. This needs testing. > There is already a field telling if you could land of not ("AMERISS") that was found by some old ARCgis coding (ARCinfo and aml, not sure exactly how). It was a consultant that did this code for us a couple of years ago, now we want to transfer to R, explaining why I have this problem now. The code was also using iteration, however not random, on rasterized lakes. I think it would be interesting to compare the result from the 2 techniques. > > I'm going to try that this afternoon, and let you know the result. I had a brief think about how to do this with rasters, but I couldn't see any advantage, unless there's some neat trick way of eroding parts of the lake that can't possibly be part of a given rectangle... Barry _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Barry Rowlingson skreiv:
>> There is already a field telling if you could land of not ("AMERISS") >> that was found by some old ARCgis coding (ARCinfo and aml, not sure >> exactly how). It was a consultant that did this code for us a couple of >> years ago, now we want to transfer to R, explaining why I have this >> problem now. The code was also using iteration, however not random, on >> rasterized lakes. I think it would be interesting to compare the result >> from the 2 techniques. >> >> I'm going to try that this afternoon, and let you know the result. > > I had a brief think about how to do this with rasters, but I couldn't > see any advantage, unless there's some neat trick way of eroding parts > of the lake that can't possibly be part of a given rectangle... Perhaps using a 2D Boyer–Moore-style algorithm on the raster could work? http://techreports.lib.berkeley.edu/accessPages/CSD-93-784.html Of course, this doesn’t help with the orientation problem, so you’ll have to repeat the algorithm (if it doesn’t succeed one the first trial) for various orientations. But, depending on the required accuracy, it might suffice with just a *few* rotations of the island or rectangle (5? 10? 20?). Or at least it should help to successively try orientations that are very different from each other (e.g., first 0° and 90°, then 45°, then 22.5° and 67.5°, etc.). -- Karl Ove Hufthammer E-mail: [hidden email] Jabber: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Karl Ove Hufthammer skreiv:
> Perhaps using a 2D Boyer–Moore-style algorithm on the raster could work? > http://techreports.lib.berkeley.edu/accessPages/CSD-93-784.html > > Of course, this doesn’t help with the orientation problem, But this one does (I think): Optimal Exact and Fast Approximate Two Dimensional Pattern Matching Allowing Rotations http://www.cwr.cl/publicaciones/cpm02.1.pdf -- Karl Ove Hufthammer E-mail: [hidden email] Jabber: [hidden email] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
In reply to this post by Karl Ove Hufthammer
On Tue, Jul 31, 2012 at 9:00 AM, Karl Ove Hufthammer <[hidden email]> wrote:
> Karl Ove Hufthammer skreiv: > >> Perhaps using a 2D Boyer–Moore-style algorithm on the raster could work? >> http://techreports.lib.berkeley.edu/accessPages/CSD-93-784.html >> >> Of course, this doesn’t help with the orientation problem, > > But this one does (I think): > > Optimal Exact and Fast Approximate Two Dimensional > Pattern Matching Allowing Rotations > http://www.cwr.cl/publicaciones/cpm02.1.pdf Typical CompSci paper - 12 pages of theorems and algebra and no actual working code! Sometimes I think CompSci people don't like the purity of their algorithms being sullied by actual implementations... Barry _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
On 31/07/12 21:36, Barry Rowlingson wrote: <SNIP> > > Typical CompSci paper - 12 pages of theorems and algebra and no > actual working code! Sometimes I think CompSci people don't like the > purity of their algorithms being sullied by actual implementations... Fortune nomination! cheers, Rolf _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
In reply to this post by Barry Rowlingson
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@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@data$NOUNIQUE[1666], lacs, iter=30000) # for one lake library(parallel) cl=makeCluster(7) system.time(is.ameri <- parSapply(cl,lacs@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 _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
In reply to this post by Barry Rowlingson
On Tue, Jul 31, 2012 at 2:21 PM,
[hidden email] <Bastien.Ferland- > 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. Considering the application, false positives are a problem! I'd rather miss out on flying to a lake that I could land on than fly to a lake that I can't! My algorithm shouldn't generate any false positives. Barry -- blog: http://geospaced.blogspot.com/ web: http://www.maths.lancs.ac.uk/~rowlings web: http://www.rowlingson.com/ twitter: http://twitter.com/geospacedman pics: http://www.flickr.com/photos/spacedman _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Free forum by Nabble | Edit this page |