inst/scripts/99.idw.gam.maps.r

In jae0/ecmgis: Package that interacts with stmv.

# p = bio.snowcrab::load.environment(p=p)

# ----------------------------
# 1 - define plot parameters

p = aegis::aegis_parameters(DS="groundfish")

set = aegis::groundfish.db( p=p, DS="set" )  ### TODO -- convert to aegis_db
set$sa = 1  # dummy required for mapping
season = "summer"
set = set[ filter.season( set$julian, period=season, index=T ) , ]

tp = lonlat2planar( set, proj.type=p$internal.crs )

grid = spatial_grid(p, DS="lonlat.coords")
tp$lon = grid_internal( tp$lon, grid$lon )
tp$lat = grid_internal( tp$lat, grid$lat )

pgrid = spatial_grid(p, DS="planar.coords")
tp$plon = grid_internal( tp$plon, pgrid$plon )
tp$plat = grid_internal( tp$plat, pgrid$plat )

tp = tp[which(tp$strat %in% 440:495),]
tp = tp[which(is.finite(tp$plon) & is.finite(tp$plat)),]


Z = bathymetry.db( p=p, DS="baseline", varnames=c("plon", "plat", "z") )  # SS to a depth of 500 m  the default used for all planar SS grids
Z = Z[which(Z$z>30),]

fn = project.datadirectory( "aegis", "data", "Science", "scotia.fundy.dat" )
a = read.table( fn, header=F)
names(a) = c('lon','lat')
a = lonlat2planar(a,proj.type='utm20')
ii = which(point.in.polygon(Z$plon,Z$plat,a$plon,a$plat) !=0)
Z = Z[ii,]

tp$depth = tp$z
tp$z = log(tp$totwgt.cod)
tp$z[ which(is.na(tp$z))] = 0
tp1 <- tp[which(tp$yr<1985),]

geog <- gstat(id='z',formula=z~1,locations= ~plon+plat,data=tp1, nmax=20,set=list(idp=1.25))
z = predict(geog,newdata=Z)


datarange = seq( 0, max(tp1$z), length.out=150)
cols = color.code( "seis", datarange )

outfn='idwcod<1985'
projectdir = tempdir()

dir.create (projectdir, showWarnings=FALSE, recursive =TRUE)
fn = file.path( projectdir, paste(outfn, "png", sep="." ) )
png( filename=fn, width=3072, height=2304, pointsize=40, res=300 )
lp = aegis_map( xyz=z[,c("plon", "plat", "z.pred")], depthcontours=T, pts=tp1[,c('plon','plat')],
            at=datarange, col.regions=cols, annot.cex=6 )
print(lp)
dev.off()




geogam <- gam(z~s(depth,k=4)+ s(plon,plat,k=200,bs='ts'),data=tp1)
names(Z)[3] = 'depth'
Z$depth = log(Z$depth)
zgam = predict(geogam,newdata=Z,type='response',se.fit=T)
zgam = cbind(Z,zgam)

datarange = seq( 0, quantile(tp1$z,0.999), length.out=150)
cols = color.code( "seis", datarange )

outfn='gamdepthcod<1985'
projectdir = tempdir()

dir.create (projectdir, showWarnings=FALSE, recursive =TRUE)
fn = file.path( projectdir, paste(outfn, "png", sep="." ) )
png( filename=fn, width=3072, height=2304, pointsize=40, res=300 )
lp = aegis_map( xyz=zgam[,c("plon", "plat", "fit")], depthcontours=T, pts=tp1[,c('plon','plat')],
            at=datarange, col.regions=cols,annot.cex=6 )
print(lp)
dev.off()


datarange = seq( 0, quantile(zgam$se.fit,0.999), length.out=150)
cols = color.code( "seis", datarange )



outfn='gamdepthcod<1985se.fit'
projectdir = tempdir()

dir.create (projectdir, showWarnings=FALSE, recursive =TRUE)
fn = file.path( projectdir, paste(outfn, "png", sep="." ) )
png( filename=fn, width=3072, height=2304, pointsize=40, res=300 )
lp = aegis_map( xyz=zgam[,c("plon", "plat", "se.fit")], depthcontours=T, pts=tp1[,c('plon','plat')],
            at=datarange, col.regions=cols,annot.cex=6 )
print(lp)
dev.off()