---
title: "Madagascar meteorological data"
author: "Duncan Golicher"
date: "13/02/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,message=FALSE,warning=FALSE)
library(giscourse)
library(tidyverse)
library(mapview)
library(dygraphs)
library(xts)
library(sf)
library(lubridate)
conn<-sconnect11("clima2")
# query<-"select id, geometry from ghcnd where st_Y(geometry)<-12 and  st_Y(geometry)>-26 and st_X(geometry)<49 and st_X(geometry)>44 "
# d<-st_read(conn,query=query)
# st_crs(d)
# mapview(d)
# f<-function(x)data.frame(id=x,get_ghcnd(x))
#   climate<-lapply(unique(d$id), f)
#   climate<-as.data.frame(do.call("rbind",climate))
# save(d,climate,file="madagascar_climate.rda")
  load("madagascar_climate.rda")
  climate %>% filter(element=="PRCP") %>% group_by(id) %>% summarise(n=n()/365) ->df
  
```

# Introduction

According to the UN Madagascar is suffering the worst drought in 40 years. 

https://news.un.org/en/story/2021/11/1104652 

https://www.bbc.co.uk/news/world-africa-58303792

It has even been claimed that the Madagascan famine is the first in history attributable to climate rather than other more avoidable factors.

https://www.independent.co.uk/climate-change/madagascar-famine-climate-crisis-b1888058.html

However the reports are based on "living memory", not on meteorological data. I downloaded meteorological records from the NOA database to check the change in rainfall and temperatures on the island. Records of daily rainfall were padded with zeros for days in which no rain was recorded. Missing data in the temperature series was filled with values from the previous day.


# Results {.tabset}

## Position of climate stations

Three climate stations have good data. The positions are shown below.

```{r}
d %>% filter(id %in% c( "MA000067083","MA000067197","MAM00067027"))->stations
mapview(stations)
```


```{r}
climate %>% filter(id %in% c( "MA000067083","MA000067197","MAM00067027")) %>% filter(element=="PRCP") %>% mutate(value=value/10) ->rain
#rain %>% group_by(id) %>% complete(date=seq.Date(from=min(date),to= max(date), by="day",fill=0)) ->rain
rain$value[is.na(rain$value)]<-0
climate %>% filter(element=="TAVG") %>%  filter(id %in% c( "MA000067083","MA000067197","MAM00067027")) %>% mutate(value=value/10)->tavg
#tavg %>% group_by(id) %>% complete(date=seq.Date(from=min(date),to= max(date), by="day",fill=0)) ->tavg
#tavg %>% fill(value)->tavg

```

```{r}
rain2<-rain
rain2$year<-lubridate::year(rain2$date)
rain2$month<-lubridate::month(rain2$date)
rain2$date<-ymd(paste(rain2$year,6,15))
rain2 %>% group_by(id,date) %>% summarise(value=sum(value)) ->rain2
```


## Rainfall in mm per day

Taking a rolling average over 365 days is a good way of visualising droughts. The rainfall in Madagascar does tend to be concentrated in a few extreme events, so changing the rolling average using the box in the bottom left corner and refreshing the figure can help to eyeball the data.

```{r}

rain$value[is.na(rain$value)]<-0
rain %>% pivot_wider(names_from = id, values_from = value) %>% arrange(date)->rn2
rn2 %>% mutate(date=as.Date(date)) ->rn2
rn2<-rn2[,-2]
rn2 %>% complete(date=seq.Date(from=min(date),to= max(date), by="day",fill=0)) ->rn2
rn2[is.na(rn2)]<-0
rn<-xts(rn2[,2:3],order.by =rn2$date)
dygraph(rn) %>% dygraphs::dyRangeSelector(dateWindow = c("1980-01-01", "2021-02-11")) %>% dyRoller(showRoller = TRUE, rollPeriod=365)                                                                           
```


## Rainfall in mm per day Mahsjang

The station in the North West of the island has been given its own graph as the climate is rather different.

```{r}


rn<-xts(rn2[,4],order.by =rn2$date)
dygraph(rn) %>% dygraphs::dyRangeSelector(dateWindow = c("1980-01-01", "2021-02-11")) %>% dyRoller(roll=365)
```



## Average temperature

```{r}

#rain$value[is.na(rain$value)]<-0
tavg %>% pivot_wider(names_from = id, values_from = value) ->rn2
rn2 %>% complete(date=seq.Date(from=min(date),to= max(date), by="day",fill=0)) ->rn2
rn2 %>% fill(MA000067083) %>% fill(MA000067197) %>% fill(MAM00067027) ->rn2
rn<-xts(rn2[,3:4],order.by =rn2$date)
dygraph(rn) %>% dygraphs::dyRangeSelector(dateWindow = c("1980-01-01", "2021-02-11")) %>% dyRoller(roll=365)
```

## Average temperature Mahsjang

```{r}

#rain$value[is.na(rain$value)]<-0
tavg %>% pivot_wider(names_from = id, values_from = value) ->rn2
rn<-xts(rn2[,5],order.by =rn2$date)
dygraph(rn) %>% dygraphs::dyRangeSelector(dateWindow = c("1980-01-01", "2021-02-11")) %>% dyRoller(roll=365)
```

## Discussion

The data does not appear to show any dramatic decrease in rainfall. In fact if anything there is a suggestion that rainfall has increased. This may not be beneficial for agriculture if the increase is the result of extreme rainfall events. However there is no clear indication from this of severe drought in any of the areas based on a consideration of the amount of rain alone. Neither is there any indication of extreme changes in temperature. Thus claiming that climate change alone is responsible for the famine requires further evidence, as it is not well supported by easily available data.


# Cyclones in Madagascar

The recent passing of cyclone Batsirai over the island has also been cited as evidence that climate change alone is responsible for the dire situation.

NOAA hold an archive of all the recorded hurricanes, typhoons, tropical cyclones and major tropical storms here.

https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/csv/

The data available on 12 February 2022 have been downloaded and imported into a PostGIS database.


```{r}

# d1 <- read_csv("ibtracs.ALL.list.v04r00.csv")[,2:12]
# nms<-tolower(names(d1))
# d1 <- read_csv("ibtracs.ALL.list.v04r00.csv",skip=1)[,2:12]
# names(d1)<-nms
# library(dplyr)
# d<- sf::st_as_sf(d1, coords = c("lon","lat"),crs=4326)
# st_write(d,con,"ibtracs")
# library(RPostgreSQL)
# RPostgreSQL::dbSendQuery(con, "CREATE INDEX idx_basin 
# ON ibtracs(basin)")
# RPostgreSQL::dbSendQuery(con, "CREATE INDEX idx_year 
# ON ibtracs(season)")
# 
```

# Cyclones since 1950 {.tabset}


```{r}
d<-st_read(conn,query="select * from ibtracs where basin like '%SI%' and season >1950")

```

## Numbers registered per year

```{r}
d %>% st_drop_geometry()  %>% group_by(season) %>% summarise(n=length(unique(number))) ->d4
ggplot(d4,aes(x=season,y=n)) +geom_line() +geom_smooth()
```

## Recorded windspeeds over 60 knots per hour

```{r}
d %>% filter(d$wmo_wind>60) -> winds
mapview(winds,zcol="wmo_wind") ->m
m@map %>% leaflet.extras::addFullscreenControl()
```

## Maximum wind speeds

```{r}
winds %>% sf::st_drop_geometry() %>% group_by(season) %>% summarise(maximum_windspeed=max(wmo_wind)) ->w2
ggplot(w2,aes(x=season,y=maximum_windspeed))+geom_line()+geom_smooth()
```

## Lowest recorded pressure

```{r}
d %>% filter(d$wmo_pres>100) -> winds
winds %>% sf::st_drop_geometry() %>% group_by(season) %>% summarise(minimum_pressure=min(wmo_pres )) ->w2
ggplot(w2,aes(x=season,y=minimum_pressure))+geom_line()+geom_smooth()
```

## Discussion

The type of storms suffered by Madagascar would not usually lead to loss of life in a developed region, but they can be deadly when people live in very precarious conditions. This is why adaptation to climate (whether it is changing or not) saves lives, livelihoods and biodiversity.
There are very few records of wind speeds over 60 knots ph over land. High winds are registered in the ocean itself, but not over most of Madagascar. The actual numbers of recorded cyclones in the Southern Indian Ocean has not been increasing recently, although some years such as 2019 have shown strong activity. Whenever looking at climate related stories it is helpful to take a long view and find all available data. Ocean warming may increase the severity of tropical depressions without increasing their numbers, so this alone does not indicate the potential for damage. An example of what could be considered to be rather misleading reporting is found here.

https://www.rmets.org/metmatters/indian-ocean-tropical-cyclones

This is a well-meaning article, which warns of the potential damage caused by cyclones and thus encourage defences to be put in place against their effects. However, the article is full of dramatic statements such as "suddenly exploded". It does not include quantitative evidence and data. The statement that "To say this is rare would be an understatement. Some meteorologists are even suggesting that this is the first time that five simultaneous tropical disturbances have ever been reported in the West Indian Ocean" could be factual, in the sense that some meteorologists could indeed suggest this. However, it is journalistic in its tone, and in its content. The statement is not clearly supported by the available data, although a deeper analysis which takes into account more details regarding the definitions of storms and measurements of wind speeds and pressure would be needed in order to clarify why it could be made. There is some evidence in the data that the recorded maximum wind speeds each season may indeed have increased since 1980. However, wind speeds over the oceans were not recorded until the satellite derived data became available. As monitoring methods improve in their temporal resolution, they become more capable of picking up the short duration peaks. So, further careful analysis is required to fully understand the data. An initial exploration does not provide strong evidence of clearly sustained recent increase in cyclone activity.
So, relying only on opinions, rather than looking at the actual data, is not advisable, even when the arguments being made are intended to enhance awareness of an important issue. The empirical data is clearly nuanced and somewhat ambiguous.

# Cyclones since 2014 {.tabset}

## Tracks as lines

```{r}

d<-st_read(conn,query="select * from ibtracs where basin like '%SI%' and season >2014")

```


```{r}

d %>% group_by(number) %>% summarise(season=min(season),number=min(number)) %>% sf::st_cast("LINESTRING") ->d3
m<-mapview(d3, burst="season", legend=FALSE) 
m@map %>% leaflet.extras::addFullscreenControl()
```

## Points with recordings

```{r}
d %>% mapview(zcol="name", legend=FALSE) ->m
m@map %>% leaflet.extras::addFullscreenControl() 

```



