0

I am using GAM to simply model a river section bathymetry from point field measurement of water depth. I've used the process several time and without problems (a great example of the workflow can be found here https://fromthebottomoftheheap.net/2016/03/27/soap-film-smoothers/). This time I am getting a wired behaviour of the smoother close to the boundary of my area, even if I forced the smooth to get a value=0 at the boundary. I see my question being similar to this old post that has no responses How to set boundary condition in a complex soap film GAM smoother?.

Here the procedure I followed:

  1. Set boundary and knots

area<-read_sf("/area.shp")
g<-st_coordinates(area)
g<-as.data.frame(g)
g$past<-paste(g$X,g$Y)

bound <- list(list(long = g[,1], lat = g[,2], f=rep(0,nrow(g))))##f= values at boundary

N <- 20#number of knots
gx <- seq(min(g[,1]), max(g[,1]), len = N)
gy <- seq(min(g[,2]), max(g[,2]), len = N)
gp <- expand.grid(gx, gy)
names(gp) <- c("long","lat")
gp_sf<-st_as_sf(gp,coords = c("long","lat"))
st_crs(gp_sf)<-4326
knots <- gp[with(gp, inSide(bound, long, lat)), ]
out<-autocruncher(bound,knots,x="long",y="lat") #download autocruncher function form https://github.com/dill/soap_checker
  1. Load filed measurement and filer out external points

bat<-read_sf("batymetry_2d.shp")[,c(3,9)] ##fied water depth points, depth in cm
bat$depth<- -1*as.numeric(bat$'_REMARKS')
bat<-bat[-which(is.na(bat$depth)|bat$depth< -100),]
bat$long<-unlist(st_geometry(bat))[seq(1,(nrow(bat)*2),2)]
bat$lat<-unlist(st_geometry(bat))[seq(2,(nrow(bat)*2),2)]

##filter out external point
bat$wt<-c(st_within(bat,area,sparse=F))
bat2<-bat[-which(bat$wt==F),]

Here the knots distribution (black points) and my measures location (blue points) enter image description here

  1. Run GAM
gam1<-gam(depth~s(long,lat,bs="so",xt=list(bnd=bound)),data=bat2
,knots = knots[-c(out),] ,method = "REML")

As you can see the model give wired estimation of water depth near the edges where they are supposed to be 0. White area inside the boundary are removal of estimation where are located of emerging stones. Measure are in cm. Depth in cm I have clearly checked any presence of outliers, and tried various number of knots but without and major improvement. I don't know if it is due to observation\knots location, but any suggestion will be appreciated.

matte_abdn
  • 31
  • 1
  • I suspect your description of what are the knots and the data is back-to-front? You are going to be doing a lot of extrapolation and forcing the data to be 0 on the boundary could well be causing the problems you are seeing. What does the spatial smooth look like if you i) don't force it to be boundary interpolating and just use a a TPRS with `s(x,y)`? and /or ii) if you don't force the boundary to be 0 but instead allow it to be estimated? Plotting your data would be useful – Gavin Simpson Nov 28 '22 at 08:03
  • Dear Gavin, thanks for this. I have added some EDITS and info. I am sorry, but I am not sure what you mean with back-to-front? thanks – matte_abdn Nov 30 '22 at 09:33
  • You say that the blue dots are your data and the block dots are the knots - hopefully that's not true? Also, you need to put the edit *into your question* not as an answer – Gavin Simpson Nov 30 '22 at 11:23
  • I asked you to plot your data - plot the depth at the points that represent where you sampled depth. Anyway, I think we've found the answer - you have to have the spatial covariates in a common format - WGS84 is like lat/long and those are not common units (I think). You need the covariates in a local coordinate system where 1 unit in the x direction is the same as 1 unit in the y direction. I hadn't noticed that you had lat/lon in your original question. – Gavin Simpson Nov 30 '22 at 11:26
  • I will also add that to use the soap film optimally, you should be including the internal boundaries around the large stones, otherwise you are smoothing through the stones, which may or may not be right... – Gavin Simpson Nov 30 '22 at 11:27
  • Thanks Gavin for all your suggestion, yes the blue dots are the knots, sorry. About the internal boundaries you are right, but at the moment I am trying to understand how to specify them without getting errors. – matte_abdn Dec 01 '22 at 09:05

0 Answers0