Input 3: The TC model spatial domain
finally, the domain and geometry for the model output needs to be
defined. The domain size and coordinates are calculated with the
land_geometry function. A domain can simply be defined with
terra::rast. Further to this a coastline polygon can be
rasterize’d to define land, and the inland distance can be
calculated with the terra::costDistance function to reduce
winds overland due to terrestrial roughness (under development and
commented out for now).
r = rast(xmin = 145,xmax=149,ymin = -19,ymax = -16.5,resolution=.01)
values(r) = 0
GEO_land = land_geometry(r,r)
#
#land_v <- vect(system.file("extdata/OSM_500m_QLD/OSM_500m_QLD.shp", package="TCHazaRds"))
#land_r = rasterize(land_v,r,touches=TRUE,background=0)
#inland_proximity = costDistance(land_r,target = 0,scale=1)
#GEO_land = land_geometry(land_r,inland_proximity)
#plot(inland_proximity,main = "Inland Distance (m)")
#plot(TC,add=TRUE)
Output Wind Feild
Now that we have the three inputs (tracks, parameters and model
output geometry) we can compute and plot the spatial wind hazard. See Making maps in R for
plotting method.
ats = seq(0, 65, length=14)
HAZi = TCHazaRdsWindField(GEO_land = GEO_land,TC = TCi,paramsTable=paramsTable)
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)
The package rasterVis:: allows pretty spatial vector
plots of the wind field via the vectorplot function (tested
on MS-Windows machine).
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"))
}
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
netcdf file for input to force hydrodynamic or wave modelling by
including outfile filename in the function call (not shown
here, see ?TCHazaRdsWindFields).
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)
Output wind time series
Time series data can be computed for a single location. Below is a
comparison of the raw IBTrACS time step and the track interpolated to 10
minute intervals.(tested on MS-Windows machine)
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)
}
Output wind Profile
Wind profiles can be calculated for a single time step. Here we
estimate the wind speed values along the profile that is 90 degrees
clockwise (at right angles) from the TC heading/bearing direction.
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)
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()
Julian O’Grady is a @csiro.au climate
scientist investigating coastal hazards and impacts.