Mapa

Load packages

if(!require("pacman")) install.packages("pacman")

## Loading required package: pacman

pacman::p_load(raster,
               sf,
               exactextractr,
               viridis,
               tidyverse,
               geobr,
               ggspatial,
               showtext,rayshader)
  
#add font
font_add_google("Open Sans", "open")
showtext_auto()

Pull municipal boundaries for Bahia

br <- read_municipality(code_muni = "BA", year = 2020, showProgress = F)

Pull gridded precipitation data

prec <- getData('worldclim', var='prec', res=10)

Calculate data frame of mean precipitation for all months

br <- cbind(br, exact_extract(prec, br, c('mean'), progress=F))

Calculate mean precipitation anual

br <- br %>% mutate(total = mean.prec1 + mean.prec2 + mean.prec3 + mean.prec4 + mean.prec5 +
                    mean.prec6 + mean.prec7 + mean.prec8 + mean.prec9 + mean.prec10 + mean.prec11 +
                    mean.prec12)

Precipitation map

p1<- ggplot(br,aes(fill = total)) +
  geom_sf(colour = alpha("gray",0.1))+
  labs(x = "Longitute", y = "Latitude", caption = "") +
  ggtitle("Bahia")+
  annotation_north_arrow(style = north_arrow_fancy_orienteering(),location="tr")+
  annotation_scale(location="br", height = unit(0.15,"cm")) +
  theme_minimal(20, base_family = "open")+
  scale_fill_viridis(option = "D", direction =  -1, "Precipitação [mm]") +
  theme(plot.title.position = "plot",plot.margin = margin(t = 25, r = 25, b = 10, l = 25))+
  guides(fill = guide_colorbar(title.position = "top", 
                               title.hjust = .5, 
                               barwidth = unit(20, "lines"), 
                               barheight = unit(.5, "lines")))+
  theme(legend.position = "bottom")+
  labs(subtitle = "Precipitação média anual", caption = "Data: WorldClim\nCid Póvoas")

p1

Precipitation by month

ba2<- crop(prec,br)
y <- mask(ba2,br)
names(y) <- month.name
plot(y, col=hcl.colors(10, rev=T))

Pull gridded altitude data

altitude <- getData('alt', country = 'BRA')
ba <- read_municipality(code_muni = "BA", year = 2020, showProgress = F)

ba <- cbind(ba, data.frame(alt=exact_extract(altitude, ba, c('majority'), progress=F)))

Altitude map of municipalities

p1<- ggplot(ba,aes(group=name_state ,fill = alt)) +
  geom_sf(colour = alpha("gray",0))+
  labs(x = "Longitute", y = "Latitude", caption = "") +
  ggtitle("Bahia")+
  annotation_north_arrow(style = north_arrow_fancy_orienteering(),location="tr")+
  annotation_scale(location="br", height = unit(0.15,"cm")) +
  theme_minimal(20, base_family = "open")+
  scale_fill_gradientn(colours = terrain.colors(10), "Elevação [m]")+
  theme(plot.title.position = "plot",plot.margin = margin(t = 25, r = 25, b = 10, l = 25))+
  guides(fill = guide_colorbar(title.position = "top", 
                               title.hjust = .5, 
                               barwidth = unit(20, "lines"), 
                               barheight = unit(.5, "lines")))+
  theme(legend.position = "bottom")+
  labs(subtitle = "Elevação", caption = "Data: GADM\nCid Póvoas")

p1

Create mask of altitude with state map

ba0 <- read_state(code_state =  "BA", year = 2020, showProgress = F)
ba1<- crop(altitude,ba0)
y <- mask(ba1,ba0)
test_spdf <- as(y, "SpatialPixelsDataFrame")
test_df <- as.data.frame(test_spdf)
colnames(test_df) <- c("value", "x", "y")

Map altitude of Bahia

ggplot()+
  geom_sf(data=ba0,aes(group=name_state),colour = alpha("gray",0))+
  geom_tile(data=test_df, aes(x=x, y=y, fill=value), alpha=0.8) + 
  ggtitle("Bahia")+
  labs(x = "Longitute", y = "Latitude", caption = "") +
  annotation_north_arrow(style = north_arrow_fancy_orienteering(),location="tr")+
  annotation_scale(location="br", height = unit(0.15,"cm")) +
  theme_minimal(20, base_family = "open")+
  scale_fill_gradientn(colours = terrain.colors(10), "Elevação [m]")+
  theme(plot.title.position = "plot",plot.margin = margin(t = 25, r = 25, b = 10, l = 25))+
  guides(fill = guide_colorbar(title.position = "top", 
                               title.hjust = .5, 
                               barwidth = unit(20, "lines"), 
                               barheight = unit(.5, "lines")))+
  theme(legend.position = "bottom")+
  labs(subtitle = "Elevação", caption = "Data: GADM\nCid Póvoas")

Digital elevation model - DEM

relev = raster_to_matrix(y)

relev %>%
  sphere_shade(texture = "imhof1") %>%
  add_water(detect_water(relev), color = "imhof1") %>%
plot_map