---
title: "Introduction to TCHazaRds"
author: "Julian O'Grady"
output: rmarkdown::html_vignette
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{Introduction_to_TCHazaRds}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

R's compatibility to easily use fast Cpp code ([Rcpp](https://github.com/RcppCore/Rcpp)) and spatial processing (e.g. [terra](https://github.com/rspatial/terra)) makes it an attractive opensource environment to study [tropical cyclones](https://en.wikipedia.org/wiki/Tropical_cyclone), aka TCs, hurricanes and typhoons. This package estimates TC vortex wind and pressure fields using parametric equations originally coded up in python by [TCRM](https://github.com/GeoscienceAustralia/tcrm) and in Cuda Cpp by [TCwindgen](https://github.com/CyprienBosserelle/TCwindgen).

TC wind fields can be computed using three model inputs of the: 1) [TC Track](#input1), 2) [Model Parameters](#input2) and 3) [Model Spatial Domain](#input3). The TCHazaRds package can be used with other visualization and spatial analysis packages to analyse the impacts of TCs.

```{r}
suppressPackageStartupMessages(require(TCHazaRds))   # this package :)
suppressPackageStartupMessages(require(terra))       # spatial analysis
suppressPackageStartupMessages(require(rasterVis))   # enhanced raster visualization https://oscarperpinan.github.io/rastervis/
suppressPackageStartupMessages(require(sp))          # spatial methods and plotting
suppressPackageStartupMessages(require(knitr))       # formatted table
suppressPackageStartupMessages(require(raster))       # convert for raster plots

```

<a name="input1"></a>

## Input 1: The TC Track  

The first thing that is required to model near- and far-field TC winds is the TC track/path. The functions in TCHazaRds require that the tracks have a "shape-file" like spatial-vector format and have attributes of pressure, date/time, location and forward speed and direction. 

```{r}
TCi = vect(cbind(c(154,154),c(-26.1,-26)),"lines",crs="epsg:4283") #track line segment
TCi$PRES = 950 #central pressure in hPa
#TCi$RMW = 40 #radius of maximum winds in km
TCi$ISO_TIME = "2022-10-04 20:00:00" #"%Y-%m-%d %H:%M:%S", tz = "UTC"
TCi$LON = geom(TCi)[1,3] #longitude
TCi$LAT = geom(TCi)[1,4] #latitude
TCi$STORM_SPD = perim(TCi)/(3*3600) #speed of the forward motion of the TC m/s
TCi$thetaFm = 90-returnBearing(TCi) #direction of the heading of the TC (Cartesian, clockwise from x axis)

```

In the above code chunk a simple track segment is defined, but historical TC tracks, e.g. from IBTRACs, can provide the input into the model. A few tracks are provided with the package, below TC Yasi is read in.

```{r}
TC <- vect(system.file("extdata/YASI/YASI.shp", package="TCHazaRds"))
TC$PRES <- TC$BOM_PRES #different agencies each provide a PRES, you need to chose one. 
TC$STORM_SPD = TC$STORM_SPD/1.94 #provided as knots, convert to m/s
TC$thetaFm = 90-returnBearing(TC) #direction of the heading of the TC (Cartesian, clockwise from x axis)
TCi = TC[46]
```

<a name="input2"></a>

## Input 2: The TC model parameters  

The second thing required to run the model is a list of parameters, which are provided for the default settings with the package and shown below.

```{r}
paramsTable = read.csv(system.file("extdata/tuningParams/defult_params.csv",package = "TCHazaRds"))
knitr::kable(paramsTable,caption = "Parameter file")
```
<a name="Input3"></a>

## 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, 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)
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 = terra::costDist(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 and Wave 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](https://r.geocompx.org/adv-map.html) for plotting method. Ocean Wave parameters can be returned with `returnWaves = TRUE`

```{r, 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]")

```

The package `rasterVis::` allows pretty spatial vector plots of the wind field via the `vectorplot` function (tested on MS-Windows machine).

```{r, 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"))
}
```

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.

```{r, 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)
```

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`). 

```{r, 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)
```

## 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) 

```{r, 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)
}
```

## 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.

```{r, 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)
```

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
```{r}
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
```{r}
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).

```{r, 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)

```


Julian O'Grady is a @csiro.au climate scientist investigating coastal hazards and impacts.