Rainfall Data Visualization from .tif File

In many development research, we will need climate data or other data from .tif file. Here I give an example of extracting data from .tif data with .shp map file. 

The rainfall data is downloaded from https://www.chc.ucsb.edu/data (Climate Hazards Center InfraRed Precipitation with Station data (CHIRPS)), and the map of Ethiopia is downloaded from https://data.humdata.org/dataset/cod-ab-eth (The Humanitarian Data Exchange). 

Here is the R code: 

library(raster)

library(sf)

library(haven)

library(tidyr)

library(ggplot2)

library(dplyr)

library(rgdal)

library(terra)

###############################################################

# read the shp file here, dsn is the folder of shp file, if your shp file is 

# at the same folder with this r file, use dsn = "." , layer is the name of shp

# file without extension

ETH_shape <- readOGR(dsn = "eth_adm_csa_bofedb_2021_shp", 

                     layer = "eth_admbnda_adm3_csa_bofedb_2021")

plot(ETH_shape) # plot it to see

ETH_map <- as(ETH_shape, "sf")

tiffpath <- the_path_of_rainfall_data.tiff 

tiff_data <- raster(tiffpath)

rainfallDataETH <- crop(tiff_data, ETH_map)

# extract the data with different districts (use the mean value)

extract_ETHdata <- raster::extract(rainfallDataETH, ETH_shape, 

                                   fun=mean, df=TRUE, na.rm = TRUE)

ETH_shapefile_cropped <- crop(ETH_shape, rainfallDataETH)

nl_data_df <- data.frame(lon = coordinates(ETH_shapefile_cropped)[, 1], 

                         lat = coordinates(ETH_shapefile_cropped)[, 2],

                         nl_mean = drop_na(extract_ETHdata))

ETH_shapefile_cropped <- crop(ETH_shape, rainfallDataETH)

plotdata = fortify(ETH_shapefile_cropped)

plotdata$id = as.integer(plotdata$id)

colnames(nl_data_df)[3] <- "id"

plotdata_final = merge(x=plotdata,y=nl_data_df,all.x=TRUE,by='id')

colnames(plotdata_final)[10] <- "Precipitation"

ggplot() +

  geom_polygon(data = plotdata_final, 

               aes(x = long, y = lat.x, group = group,fill=Precipitation),

               color="black") +

  scale_fill_gradient2(name = "Precipitation",

                       high = "green4",

                       mid = "yellow",

                       low = "orange",midpoint = 80) +

  coord_equal()  +  ggtitle("ETH Precipitation 2009.12")+ theme_void()