diff --git a/TCHazaRds_ACS.Rproj b/TCHazaRds_ACS.Rproj
index 584b854..bf73ff4 100644
--- a/TCHazaRds_ACS.Rproj
+++ b/TCHazaRds_ACS.Rproj
@@ -16,4 +16,4 @@ BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
PackageCheckArgs: --as-cran
-PackageRoxygenize: rd,collate,namespace
+PackageRoxygenize: rd,collate,namespace,vignette
diff --git a/src/TCHazaRds.dll b/src/TCHazaRds.dll
index d0bee17..6969338 100644
Binary files a/src/TCHazaRds.dll and b/src/TCHazaRds.dll differ
diff --git a/vignettes/Introduction_to_TCHazaRds.R b/vignettes/Introduction_to_TCHazaRds.R
index 149c310..fe66fba 100644
--- a/vignettes/Introduction_to_TCHazaRds.R
+++ b/vignettes/Introduction_to_TCHazaRds.R
@@ -33,8 +33,7 @@ TCi = TC[46]
## -----------------------------------------------------------------------------
paramsTable = read.csv(system.file("extdata/tuningParams/defult_params.csv",package = "TCHazaRds"))
-#paramsTable$value[6:7] = c(0,0)
-paramsTable#knitr::kable(paramsTable,caption = "Parameter file")
+knitr::kable(paramsTable,caption = "Parameter file")
## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
r = rast(xmin = 145,xmax=149,ymin = -19,ymax = -16.5,resolution=.01)
@@ -51,3 +50,88 @@ GEO_land = land_geometry(land_r,inland_proximity)
#plot(TC,add=TRUE)
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+ats = seq(0, 65, length=14)
+HAZi = TCHazaRdsWindField(GEO_land = GEO_land,TC = TCi,paramsTable=paramsTable,returnWaves = TRUE)
+library(raster) # convert for raster plots
+dummy = raster::raster()
+TC_sp = list("sp.lines",as(TC,"Spatial"),col="black")
+sp::spplot(HAZi,"Sw",at=ats,sp.layout = TC_sp,main = "Surface wind speed [m/s]")
+ats = seq(0, 16, length=9)
+sp::spplot(HAZi,"Hs0",at=ats,sp.layout = TC_sp,main = "Deep water significant wave height [m]")
+
+
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+ats = seq(0, 65, length=14)
+if (.Platform$OS.type == "windows"){
+ UV = as(c(HAZi["Uw"],HAZi["Vw"]),"Raster") #need to convert back to raster
+ rasterVis::vectorplot(UV, isField='dXY', col.arrows='white', aspX=0.002,aspY=0.002,at=ats ,
+ colorkey=list(at=ats), par.settings=viridisTheme)+latticeExtra::layer(sp.lines(as(TC,"Spatial"),col="red"))
+}
+
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+HAZ = TCHazaRdsWindFields(GEO_land=GEO_land,TC=TC,paramsTable=paramsTable)
+sp::spplot(max(HAZ$Sw),at=ats,sp.layout = TC_sp)
+
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+outdate = seq(strptime(TC$ISO_TIME[1],"%Y-%m-%d %H:%M:%S",tz="UTC"),
+ strptime(rev(TC$ISO_TIME)[1],"%Y-%m-%d %H:%M:%S",tz="UTC"),
+ 3600)
+HAZI = TCHazaRdsWindFields(outdate=outdate,GEO_land=GEO_land,TC=TC,paramsTable=paramsTable)
+sp::spplot(max(HAZI$Sw),at=ats,sp.layout = TC_sp)
+
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+
+outdate = seq(strptime(TC$ISO_TIME[1],"%Y-%m-%d %H:%M:%S",tz="UTC"),
+ strptime(rev(TC$ISO_TIME)[1],"%Y-%m-%d %H:%M:%S",tz="UTC"),
+ 600)
+
+GEO_landp = data.frame(dem=0,lons = 147,lats=-18,f=-4e-4,inlandD = 0)
+HAZts = TCHazaRdsWindTimeSereies(GEO_land=GEO_landp,TC=TC,paramsTable = paramsTable)
+HAZtsi = TCHazaRdsWindTimeSereies(outdate = outdate,GEO_land=GEO_landp,TC=TC,paramsTable = paramsTable)
+
+main = paste(TCi$NAME[1],TCi$SEASON[1],"at",GEO_landp$lons,GEO_landp$lats)
+if (.Platform$OS.type == "windows"){
+ suppressWarnings(with(HAZts,plot(date,Sw,format = "%b-%d %HZ",type="l",main = main,ylab = "Wind speed [m/s]")))
+ with(HAZtsi,lines(date,Sw,col=2))
+ legend("topleft",c("6 hrly","10 min interpolated"),col = c(1,2),lty=1)
+}
+
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+TCi$thetaFm = 90-returnBearing(TCi)
+pp <- TCProfilePts(TC_line = TCi,bear=TCi$thetaFm+90,length =150,step=1)
+#extract the GEO_land
+GEO_land_v = extract(GEO_land,pp,bind=TRUE,method = "bilinear")
+HAZp = TCHazaRdsWindProfile(GEO_land_v,TCi,paramsTable)
+
+HAZie = extract(HAZi,pp,bind=TRUE)#,method = "bilinear")
+
+wcol = colorRampPalette(c("white","lightblue","blue","violet","purple"))
+#see ?terra::plot
+plot(HAZi,"Sw",levels=ats,col = wcol(13),range = range(ats),type="continuous",all_levels=TRUE)
+#plot(HAZp,add=TRUE,cex=1.2)
+plot(HAZp,"Sw",levels=ats,col = wcol(13),range = range(ats),type="continuous",border="grey")#,all_levels=TRUE)
+lines(TC)
+
+## -----------------------------------------------------------------------------
+TCi$B = 2.2
+paramsTableCB = paramsTable
+paramsTableCB$value[paramsTableCB$param == "betaModel"] = NA
+HAZpCP = TCHazaRdsWindProfile(GEO_land_v,TCi,paramsTableCB)
+
+
+## -----------------------------------------------------------------------------
+TCi$RMAX2 = 200
+paramsTableRMAX2 = paramsTable
+paramsTableRMAX2$value[paramsTableRMAX2$param == "rMax2Model"] = NA
+HAZpRMAX2 = TCHazaRdsWindProfile(GEO_land_v,TCi,paramsTableRMAX2)
+
+
+## ----out.width = '80%',fig.height=4,fig.width=6, fig.align = "center"---------
+plot(HAZp$radialdist,HAZp$Sw,type="l",xlab = "Radial distance [km]",ylab = "Wind speed [m/s]",ylim = c(0,70));grid()
+lines(HAZp$radialdist,HAZpCP$Sw,col=2)
+lines(HAZpRMAX2$radialdist,HAZpRMAX2$Sw,col=4)
+legend("topleft",c("B = MK14, RMAX2 = 150 km",paste0("B = ",TCi$B,", RMAX2 = 150 km"),paste0("B = MK14, RMAX2 = ",TCi$RMAX2," km")),lty=1,col = c(1,2,4),cex=.7)
+title("Profiles of different peakness B and outer radius RMAX2 parameters",cex.main=.9)
+
+
diff --git a/vignettes/Introduction_to_TCHazaRds.html b/vignettes/Introduction_to_TCHazaRds.html
index 0f6fdb3..438a01e 100644
--- a/vignettes/Introduction_to_TCHazaRds.html
+++ b/vignettes/Introduction_to_TCHazaRds.html
@@ -12,7 +12,7 @@
-
+
Introduction to TCHazaRds
@@ -340,7 +340,7 @@
Introduction to TCHazaRds
Julian O’Grady
-2024-08-18
+2024-09-05
@@ -394,8 +394,7 @@ Input 2: The TC model parameters
which are provided for the default settings with the package and shown
below.
paramsTable = read.csv(system.file("extdata/tuningParams/defult_params.csv",package = "TCHazaRds"))
-#paramsTable$value[6:7] = c(0,0)
-knitr::kable(paramsTable,caption = "Parameter file")
Parameter file
@@ -414,95 +413,103 @@ Input 2: The TC model parameters
| eP |
1010.0000 |
-“environmental pressure [hPa]” |
+environmental pressure [hPa] |
| rMaxModel |
1.0000 |
-“TC radius of maximum winds model
+ | TC radius of maximum winds model
c(‘AP21’=0,‘MK14’=1,‘WR04’=2,‘VW08’=3,‘TW08’=4), or NA to use input
-TC\(RMW, see function rMax_modelsR()" |
-|vMaxModel | 1.0000|"TC maximum velocity model
-c('AP21'=0,'MK14'=1,'WR04'=2,'VW08'=3,'AH77'=4),
-or NA to use input TC\)VMAX, see function vMax_modelsR()” |
+TC$RMAX, see function rMax_modelsR
+| vMaxModel |
+1.0000 |
+TC maximum velocity model
+c(‘AP21’=0,‘MK14’=1,‘WR04’=2,‘VW08’=3,‘AH77’=4), or NA to use input
+TC$VMAX, see function vMax_modelsR |
+
+
| betaModel |
1.0000 |
-“TC beta model
+ | TC beta model
c(‘AP21’=0,‘MK14’=1,‘WR04’=2,‘VW08’=3,‘AH77’=4), or NA to use input
-TC$B, see function beta_modelsR()” |
+TC$B, see function beta_modelsR
+
+
+| rMax2Model |
+1.0000 |
+TC outer radius of 17.5m/s winds model
+(‘150km’=1,‘CK22’=2), or NA to use inuput TC$RMAX2, see function
+rMax2_modelsR |
| pressureProfileModel |
0.0000 |
-“TC pressure profile
-c(‘Holland’=0,‘McConochie’=2)” |
+TC pressure profile c(‘Holland’=0,‘McConochie’=2) |
| windProfileModel |
2.0000 |
-“TC wind profile c(‘Holland’=0,‘McConochie’=2)” |
+TC wind profile c(‘Holland’=0,‘McConochie’=2) |
| windVortexModel |
2.0000 |
-“TC wind vortex model
-c(‘Kepert’=0,‘Hubbert’=1,‘McConochie’=2,‘Jelesnianski’=4)” |
+TC wind vortex model
+c(‘Kepert’=0,‘Hubbert’=1,‘McConochie’=2,‘Jelesnianski’=4) |
| g |
9.8100 |
-“acceleration due to gravity [m/s2]” |
+acceleration due to gravity [m/s2] |
| rhoa |
1.1400 |
-“air density [kg/m3]” |
+air density [kg/m3] |
| surface |
1.0000 |
-“equals one if winds are reduced from the gradient
-level to the surface, otherwise gradient winds are returned.[-]” |
+equals one if winds are reduced from the gradient level
+to the surface, otherwise gradient winds are returned.[-] |
| Decay_a1 |
0.6150 |
-“exponential surface wind decay inland constant
-[-]” |
+exponential surface wind decay inland constant [-] |
| Decay_a2 |
0.9450 |
-“exponential surface wind decay over water constant
-[-]” |
+exponential surface wind decay over water constant
+[-] |
| Decay_a3 |
0.5125 |
-“exponential surface wind decay inland exponent
-[-]” |
+exponential surface wind decay inland exponent [-] |
| Wave_a |
0.2900 |
-“O’Grady 2024 eq.1 a parameter” |
+O’Grady 2024 eq.1 a parameter |
| Wave_x |
1.0600 |
-“O’Grady 2024 eq.1 x parameter” |
+O’Grady 2024 eq.1 x parameter |
| Wave_b |
-0.0157 |
-“O’Grady 2024 eq.1 b parameter” |
+O’Grady 2024 eq.1 b parameter |
| Wave_c |
-0.0294 |
-“O’Grady 2024 eq.1 c parameter” |
+O’Grady 2024 eq.1 c parameter |
@@ -543,10 +550,10 @@ Output Wind and Wave Feild
dummy = raster::raster()
TC_sp = list("sp.lines",as(TC,"Spatial"),col="black")
sp::spplot(HAZi,"Sw",at=ats,sp.layout = TC_sp,main = "Surface wind speed [m/s]")
-
+
ats = seq(0, 16, length=9)
sp::spplot(HAZi,"Hs0",at=ats,sp.layout = TC_sp,main = "Deep water significant wave height [m]")
+
The package rasterVis:: allows pretty spatial vector
plots of the wind field via the vectorplot function (tested
on MS-Windows machine).
@@ -556,14 +563,14 @@ Output Wind and Wave Feild
rasterVis::vectorplot(UV, isField='dXY', col.arrows='white', aspX=0.002,aspY=0.002,at=ats ,
colorkey=list(at=ats), par.settings=viridisTheme)+latticeExtra::layer(sp.lines(as(TC,"Spatial"),col="red"))
}
-
+
The hazard can be also calculated for the entire track too (by adding
a s to the end of TCHazaRdsWindField to make
it plural), and then the maximum wind speed at each grid cell can be
plotted.
HAZ = TCHazaRdsWindFields(GEO_land=GEO_land,TC=TC,paramsTable=paramsTable)
sp::spplot(max(HAZ$Sw),at=ats,sp.layout = TC_sp)
+
The track can be interpolate to say, hourly intervals by defining an
outdate from the start to the end date of the TC, stepping
by 3600 seconds. The output from these functions can be written to a
@@ -575,7 +582,7 @@
Output Wind and Wave Feild
3600)
HAZI = TCHazaRdsWindFields(outdate=outdate,GEO_land=GEO_land,TC=TC,paramsTable=paramsTable)
sp::spplot(max(HAZI$Sw),at=ats,sp.layout = TC_sp)
-
+
Output wind time series
@@ -597,7 +604,7 @@
Output wind time series
with(HAZtsi,lines(date,Sw,col=2))
legend("topleft",c("6 hrly","10 min interpolated"),col = c(1,2),lty=1)
}
-
+
Output wind Profile
@@ -615,15 +622,33 @@
Output wind Profile
wcol = colorRampPalette(c("white","lightblue","blue","violet","purple"))
#see ?terra::plot
plot(HAZi,"Sw",levels=ats,col = wcol(13),range = range(ats),type="continuous",all_levels=TRUE)
-
+
#plot(HAZp,add=TRUE,cex=1.2)
plot(HAZp,"Sw",levels=ats,col = wcol(13),range = range(ats),type="continuous",border="grey")#,all_levels=TRUE)
lines(TC)
+
+TC wind fields can be modelled, or tested, with observed, or
+constant, B (Beta) profile peakedness parameter by defining TC$B and
+setting betaModel = NA in paramsTable
+TCi$B = 2.2
+paramsTableCB = paramsTable
+paramsTableCB$value[paramsTableCB$param == "betaModel"] = NA
+HAZpCP = TCHazaRdsWindProfile(GEO_land_v,TCi,paramsTableCB)
Other parameters can be adjusted, here we model a larger outer radius
+(RMAX2) profile parameter by defining TC$RMAX2 and setting rMax2Model =
+NA in paramsTable
+TCi$RMAX2 = 200
+paramsTableRMAX2 = paramsTable
+paramsTableRMAX2$value[paramsTableRMAX2$param == "rMax2Model"] = NA
+HAZpRMAX2 = TCHazaRdsWindProfile(GEO_land_v,TCi,paramsTableRMAX2)
Positive radial distance values are to the right of the forward
motion (90 deg clockwise).
-plot(HAZp$radialdist,HAZp$Sw,type="l",xlab = "Radial distance [km]",ylab = "Wind speed [m/s]");grid()
+plot(HAZp$radialdist,HAZp$Sw,type="l",xlab = "Radial distance [km]",ylab = "Wind speed [m/s]",ylim = c(0,70));grid()
+lines(HAZp$radialdist,HAZpCP$Sw,col=2)
+lines(HAZpRMAX2$radialdist,HAZpRMAX2$Sw,col=4)
+legend("topleft",c("B = MK14, RMAX2 = 150 km",paste0("B = ",TCi$B,", RMAX2 = 150 km"),paste0("B = MK14, RMAX2 = ",TCi$RMAX2," km")),lty=1,col = c(1,2,4),cex=.7)
+title("Profiles of different peakness B and outer radius RMAX2 parameters",cex.main=.9)
Julian O’Grady is a @csiro.au climate
scientist investigating coastal hazards and impacts.