Skip to content

dill/soap_checker

master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Code

Latest commit

 

Git stats

Files

Permalink
Failed to load latest commit information.
Type
Name
Latest commit message
Commit time
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Check whether a soap film smoother boundary and knots make sense

Setting up a soap film smoother is often a hard and frustrating process. This function checks that the data that you feed to soap are “correct”. Feed it our boundary and optionally knots and data and the function will show you what soap will model.

The function also plots the area that will be modelled by the soap film smoother in red.

To do all the spatial stuff, we require a couple of extra libraries: rgeos and sp (as well as mgcv).

Ensuring the boundary is right

Simple example with the “Ramsay horseshoe”

source("soap_check.R")
fsb <- list(fs.boundary())
soap_check(fsb)

## [1] TRUE

“inside-out” "Ramsay horseshoe

fsb_io <- fs.boundary()
fsb_io <- list(fsb_io,
               list(x=range(fsb_io$x)[c(1,1,2,2,1)],
                    y=range(fsb_io$y)[c(1,2,2,1,1)]))
soap_check(fsb_io)

## [1] TRUE

Knots

Sometimes knots are too close to the boundary, this end up with errors like:

Error in crunch.knots(ret$G, knots, x0, y0, dx, dy) :
  knot 10 is on or outside boundary

Generated by this code:

## Simulate some fitting data, inside boundary...
set.seed(0)
n<-600
v <- runif(n)*5-1;w<-runif(n)*2-1
y <- fs.test(v,w,b=1)
names(fsb[[1]]) <- c("v","w")
ind <- inSide(fsb,x=v,y=w) ## remove outsiders
y <- y + rnorm(n)*.3 ## add noise
dat <- data.frame(y = y[ind],
                  v = v[ind],
                  w = w[ind])
# setup knots
knots <- expand.grid(v=seq(min(fsb[[1]]$v), max(fsb[[1]]$v), len=15),
                     w=seq(min(fsb[[1]]$w)+0.05, max(fsb[[1]]$w), len=10))
v <- knots$v
w <- knots$w
ind <- inSide(fsb,x=v,y=w) ## remove outsiders
knots <- knots[ind, ]
b <- gam(y~s(v,w,k=30,bs="so",xt=list(bnd=fsb)),knots=knots, data=dat)

One can go through, one-by-one and remove the offending knots as crunch.knots finds them, but that’s a bit tedious. Enter autocruncher, which returns all of the offending knots at once:

source("autocrunch.R")
crunch_ind <- autocruncher(fsb, knots, k=30, xname="v", yname="w")
crunch_ind
## [1]  1 13 66 93

We can see these knots are on the boundary are an issue (crosses):

plot(fsb[[1]]$v, fsb[[1]]$w, type="l", asp=1, xlab="v", ylab="w")
points(knots)
points(knots[crunch_ind, ], pch=4)

We can now simply remove them and fit our model successfully (via knots[-crunch_ind, ]).

Note that you need to set the k argument to autocruncher to be the same as your planned value in gam and if you set the nmax parameter for the soap film you need to provide that to autocruncher too.

Other tips

  • The script assumes that your spatial variable names are x and y.
  • Make sure that your locations are in Northings/Eastings. Using latitude and longitude will give strange results (as the soap film smoother is isotropic so treats 1 unit change in either dimension is equal, this isn’t true for lat/long!).
  • Sometimes you need to increase the tolerance (e.g tol=1e-6)
  • Note that the boundary must be a list of lists or data.frames. So if you have a polygon boundary with your boundary vertices in it, that must be wrapped in a list!

Written by David L Miller and released under the GPL (version 2).

About

Test that a soap film smoother boundary makes sense before running the model...

Resources

License

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published

Languages