require(terra)
dem <- rast(system.file("extdata/DEMs/YASI_dem.tif", package="TCHazaRds"))
land <- dem; land[land > 0] = 0
inland_proximity = distance(land,target = 0)
GEO_land = land_geometry(dem,inland_proximity)
TCi = vect(cbind(c(154,154),c(-26.1,-26)),"lines",crs="epsg:4283") #track line segment
TCi$PRES = 950
TCi$RMAX = 40
TCi$VMAX = 60
TCi$B = 1.4
TCi$ISO_TIME = "2022-10-04 20:00:00"
TCi$LON = geom(TCi)[1,3]
TCi$LAT = geom(TCi)[1,4]
TCi$STORM_SPD = perim(TCi)/(3*3600) #m/s
TCi$thetaFm = 90-returnBearing(TCi)
#OR
TC <- vect(system.file("extdata/YASI/YASI.shp", package="TCHazaRds"))
TC$PRES <- TC$BOM_PRES
TCi = TC[47]
TCi$thetaFm = 90-returnBearing(TCi)
#extract a profile/transect at right angles (90 degrees) from the TC heading/bearing direction
pp <- TCProfilePts(TC_line = TCi,bear=TCi$thetaFm+90,length =100,step=1)
#plot(dem);lines(TCi,lwd = 4,col=2)
#points(pp)
GEO_land_v = extract(GEO_land,pp,bind=TRUE,method = "bilinear")
paramsTable = read.csv(system.file("extdata/tuningParams/defult_params.csv",package = "TCHazaRds"))
#calculate the wind hazard
HAZ = TCHazaRdsWindProfile(GEO_land_v,TCi,paramsTable)
#plot(HAZ$radialdist,HAZ$Sw,type="l",xlab = "Radial distance [km]",ylab = "Wind speed [m/s]");grid()
#plot(HAZ,"Sw",type="continuous")
Run the code above in your browser using DataLab