determining suitable lakes for landing

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

determining suitable lakes for landing

Bastien.Ferland-Raymond
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
Reply | Threaded
Open this post in threaded view
|

Re: determining suitable lakes for landing

Barry Rowlingson
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
Reply | Threaded
Open this post in threaded view
|

Re: determining suitable lakes for landing

Barry Rowlingson
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
Reply | Threaded
Open this post in threaded view
|

RE : determining suitable lakes for landing

Bastien.Ferland-Raymond
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
Reply | Threaded
Open this post in threaded view
|

Re: RE : determining suitable lakes for landing

Barry Rowlingson
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
Reply | Threaded
Open this post in threaded view
|

Re: RE : determining suitable lakes for landing

Karl Ove Hufthammer
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
Reply | Threaded
Open this post in threaded view
|

Re: RE : determining suitable lakes for landing

Karl Ove Hufthammer
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
Reply | Threaded
Open this post in threaded view
|

Re: RE : determining suitable lakes for landing

Barry Rowlingson
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
Reply | Threaded
Open this post in threaded view
|

Re: RE : determining suitable lakes for landing

Rolf Turner

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
Reply | Threaded
Open this post in threaded view
|

RE : RE : determining suitable lakes for landing

Bastien.Ferland-Raymond
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
Reply | Threaded
Open this post in threaded view
|

Re: RE : RE : determining suitable lakes for landing

Barry Rowlingson
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