rm(list=ls())
setwd("C:/Dropbox/Research/_DeLoach/Ag_Data_News/")

#############################################################
## Code to create plots in Ag Data News article            ##
##  Title:          ##
##  Link:  https://asmith.ucdavis.edu/news/    ##
#############################################################


# If a package is installed, it will be loaded. If any are not, the missing package(s) will be installed 
# from CRAN and then loaded.
pacman::p_load(ggplot2,readxl,dplyr,patchwork,hrbrthemes,httr,RColorBrewer,ggpubr,tidyr,ggtext)


prod <- read.csv("https://files.asmith.ucdavis.edu/FAOSTAT_production.csv") %>%
  mutate(Country=replace(Country, Country=="United States of America", "USA"))

############################################
# Top Production
top_5_countries <- prod %>%
  filter(Element=="Production") %>%
  group_by(Country) %>%
  summarise(Production_mean=mean(Value, na.rm=T)) %>%
  arrange(-Production_mean) %>%
  slice(1:5) %>%
  pull(Country)
  
prod_propotion <- prod %>%
  filter(Element=="Production") %>%
  group_by(Country,Year) %>%
  summarise(Production_mean=mean(Value, na.rm=T)) %>%
  mutate(Country=ifelse(Country %in% top_5_countries, Country, "Other")) %>%
  group_by(Country,Year) %>%
  summarise(Production_mean=sum(Production_mean, na.rm=T))%>%
  ggplot(aes(x=Year, y=Production_mean/1000000, fill=factor(Country,levels=c(top_5_countries,"Other"))))+ 
    theme_minimal()+
    #geom_bar(width = 1, stat = "identity")+
    geom_area()+
    labs(fill="Country",y="Million Tonne", x="", title="Global Average Annual Rice Production \n (2000-2020)", caption = "https://agdatanews.substack.com \n Source: FAOSTAT, https://www.fao.org/faostat/en/#data/QCL")+ 
    theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 1, size = 8), text = element_text(size=14))
prod_propotion
ggsave("prod_propotion.png",bg="white")


######################################
# price
inflaction <- read.csv("https://files.asmith.ucdavis.edu/WorldBank_Inflaction_rate.csv")
india_price <- read.csv("https://files.asmith.ucdavis.edu/Rice_Prices_India_INR.csv") %>%
  group_by(Year) %>%
  summarise(price=mean(priceperkg))
india_price <-  merge(india_price,inflaction, by="Year") 
india_price <-  india_price %>% mutate(india.real.price = price/India.Index*100)
price <- read.csv("https://files.asmith.ucdavis.edu/FAO_PI.csv")
usa_price <-read.csv("https://files.asmith.ucdavis.edu/FAOSTAT_price.csv") %>%
  filter(Area=="United States of America") %>%
  select(Year, Value) 
usa_price <- merge(usa_price, inflaction, by="Year")
usa_price <- usa_price %>% mutate(usa.real.price = Value/USA.Index*100)
taiwan <- read.csv("https://files.asmith.ucdavis.edu/Taiwan.csv")
caption2 <- "The base period is 2003 and is adjusted into real terms by inflation rates in the countries and the world average, respectively. Source: FAO (world and USA prices), DGCIS (Indian prices). Agricultural statistics yearbook of Taiwan (Taiwan prices), World Bank (Inflation rates in India, USA, and the world), DGBAS (Inflation rates in Taiwan). https://agdatanews.substack.com"
merge_price <- merge(price, taiwan, by="Year")
merge_price <- merge(merge_price, usa_price, by="Year")
merge_price <- merge(merge_price, india_price, by="Year") %>%
  mutate(USA=usa.real.price/234.5464*100,
         Taiwan=taiwan.real.price/18.88301*100,
         World=Real.Index/77.86692*100,
         India=india.real.price/28.07293*100) %>%
  select("Year","USA","Taiwan","India","World") %>%
  pivot_longer(!Year, names_to = "Country", values_to = "PriceIndex")

price_rate <- merge_price %>% ggplot() + 
  geom_line(aes(x=Year, y=PriceIndex, color=Country)) +
  theme_minimal()+
  labs(y="Price Index", title="Rice Producer and World Export Price Index in Real Terms \n (2003=100)", caption = str_wrap(caption2, 140)) +
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 0, size = 8), text = element_text(size=14))
price_rate
ggsave("price_rate.png",bg="white")

#################################################################
# Area rate
caption1 <- "The change rate is defined as the percentage different from last year. \n Source: FAOSTAT, https://www.fao.org/faostat/en/#data/QCL. \n https://agdatanews.substack.com "
area_rate <- prod %>%
  filter(Element=="Area harvested" & Country %in% c("India", "Taiwan", "USA")) %>%
  mutate(rate=100*(Value-lag(Value))/lag(Value)) %>%
  filter(Year!="2000") %>%
  ggplot(aes(x=Year, y=rate, colour=Country))+
    theme_minimal()+
    labs(y="Percent Change", title="Percent Change in Annual Rice Growing Area \n (2001-2020)", caption = caption1)+ 
    theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 1, size = 8), text = element_text(size=14))+
    geom_line()
area_rate
ggsave("area_rate.png",bg="white")

##################################################
# Indian production
india_prod <- read.csv("https://files.asmith.ucdavis.edu/FAO_India_Rice.csv") %>%
  filter(Element == "Production", Year>=2000, Item =="Rice") %>%
  ggplot() +
  geom_line(aes(x=Year, y=Value/1000000))+
  theme_minimal()+
  labs(y="Million Tonne", x="")+ 
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 1, size = 8))
india_area <- read.csv("https://files.asmith.ucdavis.edu/FAO_India_Rice.csv") %>%
  filter(Element == "Area harvested", Year>=2000, Item =="Rice") %>%
  ggplot() +
  geom_line(aes(x=Year, y=Value/1000))+
  theme_minimal()+
  labs(y="Thousand Ha", x="", title = "Rice Production in India \n (2000-2020) \n \n")+ 
  ylim(30000,50000)+
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 1, size = 8))
india_yield <- read.csv("https://files.asmith.ucdavis.edu/FAO_India_Rice.csv") %>%
  filter(Element == "Yield", Year>=2000, Item =="Rice") %>%
  ggplot() +
  geom_line(aes(x=Year, y=Value/10))+
  theme_minimal()+
  labs(y="Kg/Ha", x="", caption = "Source: FAOSTAT, https://www.fao.org/faostat/en/#data/QCL \n https://agdatanews.substack.com")+ 
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 1, size = 8))

india <- ggarrange(india_prod, india_area, india_yield, 
                   labels=c("Production", "Area Harvested", "Yield"),
                   align = "hv",
                   vjust = 6.5,
                   font.label= list(face="bold", size=11),
                   ncol=3)
india
ggsave("india.png",bg="white")


################################################
# export
india_export <- read.csv("https://files.asmith.ucdavis.edu/FAO_India_Trade.csv") %>%
  filter(Item %in% c("Rice, milled", "Rice, broken"), Year>=2000,Element=="Export Quantity")%>%
  ggplot()+
  geom_line(aes(y=Value/1000, x=Year, color=Item))+
  theme_minimal()+
  labs(y="Thousand Tonne", x="", title = "Rice Exports from India \n (2000-2020)",caption = "Source: FAOSTAT, https://www.fao.org/faostat/en/#data/TCL  \n https://agdatanews.substack.com")+ 
  theme(plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 1, size = 8), text = element_text(size=14))
india_export  
ggsave("india_export.png",bg="white")