rainfall.Rmd

---
title: "Rainfall"
output:
  html_document:
    toc: true
    toc_float: true
    number_sections: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(tidyverse)
library(DT)
library(terra)
library(rnaturalearth)
library(viridis)
library(rgeoboundaries)
```

#### Data set details

|||
| ----------- | ----------- |
| **Data set description:** | Rainfall data worldwide and in China |
| **Source:**   | [NASA POWER](https://power.larc.nasa.gov/): The POWER Project |
| **Details on the retrieved data:**   | Monthly and yearly rainfall in Australia. |
| **Spatial and temporal resolution:**   | Hourly, daily, monthly, or yearly data for regions defined by the `lonlat` parameter. Each cell corresponds to 1/2 x 1/2 degree. |

## Downloading rainfall data with `nasapower`
 
The `nasapower` package aims at making it quick and easy to automate downloading [NASA POWER](https://power.larc.nasa.gov/) (NASA Prediction of Worldwide Energy Resource) global meteorology, surface solar energy and climatology data.


## Installing `nasapower` package 

The `nasapower` package can be directly downloaded from CRAN as follows.

```{r eval=FALSE}
install.packages("nasapower")
```

Loading the `nasapower` package:

```{r loading}
library(nasapower)
```

## Using `get_power()` to fetch data

The `get_power()` function has five arguments and returns a data frame with a metadata header in the current R session.
It has the following arguments:

- `pars`: character vector of solar, meteorological or climatology parameters to download.

- `community`: character vector providing community name. Supported values are `"AG"`, `"SB"` or `"RE"`.

`"AG"`: provides access to the agroclimatology archive, which contains industry-friendly parameters for input to crop models.
    
`"SB"`: provides access to the sustainable buildings archive, which contains parameters for the building community.
    
`"RE"`: provides access to the renewable energy archive, which contains parameters very specific to assist in the design of solar and wind powered renewable energy systems.
    
- `temporal_api`: temporal API end-point for data being queried. Supported values are `"HOURLY"`, `"DAILY"`, `"MONTHLY"`, or `"CLIMATOLOGY"`.
  
`"HOURLY"`: hourly average of `pars` by hour, day, month and year.

`"DAILY"`: daily average of `pars` by day, month and year.

`"MONTHLY"`: monthly average of `pars` by month and year.
    
`"CLIMATOLOGY"`: provide parameters as 22-year climatologies (solar) and 30-year climatologies (meteorology); the period climatology and monthly average, maximum, and/or minimum values.
    
- `lonlat`: numeric vector of geographic coordinates for a cell or region or `"GLOBAL"` for global coverage.

For single point supply a length-two numeric vector giving the decimal degree longitude and latitude in that order for the data to download. 
    
For regional coverage supply a length-four numeric as lower left (lon, lat) and upper right (lon, lat) coordinates as `lonlat = c(xmin, ymin, ymax, ymax)`
    
- `dates`: start and end dates. If only one date is provided, it will be treated as both the start and the end date and only a day's values will be returned. 
This argument should not be used when `temporal_api` is set to `"CLIMATOLOGY"`.

To know the different weather values from POWER provided within this function, type `?get_power`, and in the arguments section, click on the highlighted parameters, which goes to a page which has all the available parameters.

To download rainfall, we use the `pars = "PRECTOTCORR"`. 

### Fetching daily data for single point
   
```{r getdata}
data <- get_power(
  community = "RE",
  lonlat = c(134.489563, -25.734968),
  dates = c("2000-01-01", "2000-05-01"),
  temporal_api = "DAILY",
  pars = "PRECTOTCORR"
)

data %>% datatable(extensions = c("Scroller", "FixedColumns"), options = list(
  deferRender = TRUE,
  scrollY = 350,
  scrollX = 350,
  dom = "t",
  scroller = TRUE,
  fixedColumns = list(leftColumns = 3)
))
```

### Fetching daily data for an area

```{r daily}
daily_area <- get_power(
  community = "AG",
  lonlat = c(150.5, -28.5, 153.5, -25.5),
  pars = "PRECTOTCORR",
  dates = c("2004-09-19", "2004-09-29"),
  temporal_api = "DAILY"
)

daily_area %>% datatable(extensions = c("Scroller", "FixedColumns"), options = list(
  deferRender = TRUE,
  scrollY = 350,
  scrollX = 350,
  dom = "t",
  scroller = TRUE,
  fixedColumns = list(leftColumns = 3)
))
```

### Fetching climatology data

Using `lonlat` equivalent to Australia, we will set `temporal_api = "CLIMATOLOGY"`. Since the maximum processed area is 4.5 x 4.5 degrees (100 points), we will do this in step-by-step.

```{r climate}
flag <- 1
for (i in seq(110, 150, 5)) {
  for (j in seq(-40, -10, 5)) {
    climate_avg_temp <- get_power(community = "AG",
                                  pars = "PRECTOTCORR",
                                  lonlat = c(i, (j - 5), (i + 5), j),
                                  temporal_api = "CLIMATOLOGY")
    if (flag == 1) {
      climate_avg <- climate_avg_temp
      flag <- 0
    } else{
      climate_avg <- rbind(climate_avg, climate_avg_temp)
    }
  }
}

climate_avg %>% datatable(extensions = c("Scroller", "FixedColumns"), options = list(
  deferRender = TRUE,
  scrollY = 350,
  scrollX = 350,
  dom = "t",
  scroller = TRUE,
  fixedColumns = list(leftColumns = 3)
))
```

Next, we will plot the retrieved data.

## Creating a map of annual rainfall using all data retrieved

```{r}
library(rnaturalearth)
library(raster)

# Getting world map
map <- ne_countries(country = 'australia', returnclass = "sf")

# Converting data to raster
r <- rasterFromXYZ(climate_avg[, c("LON", "LAT", "ANN")])

# Converting the raster into a data.frame
r_df <- as.data.frame(r, xy = TRUE, na.rm = TRUE)

# Plot
ggplot() +
  geom_raster(data = r_df, aes(x = x, y = y, fill = ANN)) +
  geom_sf(data = map, inherit.aes = FALSE, fill = NA) +
  scale_fill_viridis() +
  labs(
    title = "Rainfall in inches",
    fill = "Annual Rainfall",
    subtitle = "Annual rainfall in Australia"
  ) + 
  labs(x = "Longitude", y = "Latitude")
```

## Creating maps of monthly rainfall

```{r}
r <- list()
for (k in colnames(climate_avg)[-c(1:3, 16)]) {
  r[[k]] <- rasterFromXYZ(climate_avg[, c("LON", "LAT", k)])
}
r <- stack(r)
plot(r, asp = 1)
```




## References

- `nasapower` package: https://github.com/ropensci/nasapower

- NASAPOWER project: https://power.larc.nasa.gov/


---

<p style ="color:grey; font-size: 13px;">  
Last updated: `r Sys.Date()`  
Source code: https://github.com/rspatialdata/rspatialdata.github.io/blob/main/rainfall.Rmd 
<details>
<summary><span style ="color:grey; font-size: 13px;">Tutorial was complied using: (click to expand)</span></summary>
```{r echo=FALSE}
sessionInfo()
```
</details>
</p>