Mapa
Load packages
if(!require("pacman")) install.packages("pacman")## Loading required package: pacmanpacman::p_load(raster,
sf,
exactextractr,
viridis,
tidyverse,
geobr,
ggspatial,
showtext,
terra,
rayshader,
RColorBrewer,
ggpubr)
#add font
font_add_google("Open Sans", "open")
showtext_auto()
# Set config
STATE <- "SP" #State code with two letter
ESTADO <- "São Paulo" #State name Download State boundaries, from IBGE
state <- read_municipality(code_muni = STATE, year = 2020, showProgress = F)Download world clim rain data from Univ. of California (Campus of Davis)
prec <- getData('worldclim', var='prec', res=2.5)Calculate data frame of mean precipitation for all months
state <- cbind(state, exact_extract(prec, state, c('mean'), progress=F))Crop rain data using State boundary and create raster
state1 <- crop(prec,state)
y <- mask(state1,state)
spdf <- as(y, "SpatialPixelsDataFrame")
spdf <- as.data.frame(spdf)Calc annual total rain data
spdf <- spdf %>% mutate(total = prec1 +
prec2 +
prec3 +
prec4 +
prec5 +
prec6 +
prec7 +
prec8 +
prec9 +
prec10 +
prec11 +
prec12) %>%
select(total,x,y)
colnames(spdf) <- c("value", "x", "y")Plot Avg annual total rain
p1<- ggplot()+
geom_sf(data=state,aes(group=name_state),colour = alpha("gray",0))+
geom_tile(data=spdf, aes(x=x, y=y, fill=value), alpha=0.8) +
ggtitle(ESTADO)+
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_distiller(palette = "Blues",direction = 1, name = "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(25, "lines"),
barheight = unit(.4, "lines")))+
theme(legend.position = "bottom")+
labs(subtitle = "Precipitação média anual", caption = "Data: WorldClim\nCid Póvoas")
p1
Calculate mean precipitation anual
state <- state %>% 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)Mean precipitation anual map by counties
p2<- ggplot(state,aes(fill = total)) +
geom_sf(colour = alpha("gray",0.1))+
labs(x = "Longitute", y = "Latitude", caption = "") +
ggtitle(ESTADO)+
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_distiller(palette = "Blues",direction = 1, name = "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 por municÃpio", caption = "Data: WorldClim\nCid Póvoas")
p2
#
Precipitation by month
state2<- crop(prec,state)
y2 <- mask(state2,state)
names(y2) <- month.name
plot(y2, col=brewer.pal(n = 9, name = "Blues"))
Pull gridded altitude data
the data were aggregated from SRTM 90m resolution
altitude <- getData('alt', country = 'BRA',mask=TRUE)
altitude2 <- altitude
altitude <- rast(altitude * 1)
state3 <- read_municipality(code_muni = STATE, year = 2020, showProgress = F)## Using year 2020state3 <- st_transform(state3, crs = 4326)
state3 <- cbind(state3, data.frame(alt=exact_extract(altitude, state3, c('majority'), progress=F)))Altitude map of counties
p4<- ggplot(state3,aes(group=name_state ,fill = alt)) +
geom_sf(colour = alpha("gray",0))+
labs(x = "Longitute", y = "Latitude", caption = "") +
ggtitle(ESTADO)+
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")
p4
Create mask of altitude with state map
state4<- crop(altitude2,state3)
y3 <- mask(state4,state3)
spdf2 <- as(y3, "SpatialPixelsDataFrame")
spdf2 <- as.data.frame(spdf2)
colnames(spdf2) <- c("value", "x", "y")Map altitude of São Paulo
p5 <- ggplot()+
geom_sf(data=state3,aes(group=name_state),colour = alpha("gray",0))+
geom_tile(data=spdf2, aes(x=x, y=y, fill=value), alpha=0.8) +
ggtitle(ESTADO)+
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")
p5 
Calculate slope
area_slope <- terrain(y3, opt = 'slope', unit = 'degrees')
spdf3 <- as(area_slope, "SpatialPixelsDataFrame")
spdf3 <- as.data.frame(spdf3)
colnames(spdf3) <- c("value", "x", "y")spdf3 %>% mutate(value=value/100) %>% ggboxplot(y="value",
xlab = "",
ylab="Slope value",
title = "Slope",
fill = "steelblue",
bxp.errorbar = T)+
scale_y_continuous(labels = scales::percent)+
coord_flip()+
theme_minimal(20, base_family = "open")+
theme(legend.text = element_text(size=15),
axis.text.y = element_blank())
spdf3$value %>% summary## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.6355 0.9893 1.4433 1.4920 25.2814Slope map
p6<- ggplot()+
geom_sf(data=state,aes(group=name_state),colour = alpha("gray",0))+
geom_tile(data=spdf3, aes(x=x, y=y, fill=value)) +
ggtitle(ESTADO)+
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 = rev(terrain.colors(16)), "Declividade [%]")+
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 = "Declividade", caption = "Data: GADM\nCid Póvoas")
p6
calculate terrain aspect
area_aspect <- terrain(y3, opt = 'aspect', unit = 'degrees')
spdf4 <- as(area_aspect, "SpatialPixelsDataFrame")
spdf4 <- as.data.frame(spdf4)
colnames(spdf4) <- c("value", "x", "y")
spdf4 <- spdf4 %>% mutate(pcat = cut(value, c(-1,
0,
22.5,
67.5,
112.5,
157.5,
202.5,
247.5,
292.5,
337.5,
360),
label=c("Flat (-1)",
"North (0°to 22.5°)",
"Northeast (22.5° to 67.5°)",
"East (67.5° to 112.5°)",
"Southeast (112.5° to 157.5°)",
"South (157.5° to 202.5°)",
"Southwest (202.5° to 247.5°)",
"West (247.5° to 292.5°)",
"Northwest (292.5° to 337.5°)",
"North (337.5° to 360°)")))Terrain aspect map
p7<- ggplot()+
geom_sf(data=state,aes(group=name_state),colour = alpha("gray",0))+
geom_tile(data=spdf4, aes(x=x, y=y, fill=pcat), alpha=0.8) +
ggtitle(ESTADO)+
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_manual(values = rainbow(10), "")+
theme(plot.title.position = "plot",plot.margin = margin(t = 25, r = 25, b = 10, l = 25))+
theme(legend.position = "bottom",
legend.text = element_text(size=15))+
labs(subtitle = "Aspecto do terreno", caption = "Data: GADM\nCid Póvoas")
p7
Digital elevation model - DEM
relev = raster_to_matrix(y3)
relev %>%
sphere_shade(texture = "imhof1") %>%
add_water(detect_water(relev), color = "white") %>%
plot_map