1. Introduction

library (tidyverse)
library(dplyr)
library(ggplot2)
library(sf)

3. Readinig a table with Santander´s agriculture statucs

eva_santander <- read_csv(file="C:/users/morea/Documents/UNAL/2021 - 1S/GeomaticaBasica/eva_santander_c.csv")

4. Reading a shapefile with Santander´s municipalities

mun_santander<-st_read("C:/users/morea/Documents/UNAL/2021 - 1S/GeomaticaBasica/MGN2017_68_SANTANDER 2/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\morea\Documents\UNAL\2021 - 1S\GeomaticaBasica\MGN2017_68_SANTANDER 2\68_SANTANDER\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## geographic CRS: WGS 84
mun_santander
eva_santander

Let´s calculate the area of each municipality in plain coordinates (Km^2)

mun_santander$KM2 <- st_area(st_transform(mun_santander,3116))/1E6 #Lo divido por 10^6 y me transforma de m a km2
mun_santander$KM2 <-as.numeric(mun_santander$KM2)
mun_santander$KM2 <- round(mun_santander$KM2,3)
min(mun_santander$KM2)
## [1] 19.694
max(mun_santander$KM2)
## [1] 3174.289

Let’s make a map of municipalities:

library(leaflet)
bins <- c(0, 150, 300, 450, 600, 750, 900, 1200, 1600, 2000, 3200)
pal <- colorBin("RdYlGn", domain = mun_santander$KM2, bins = bins)

  mapa <- leaflet(data = mun_santander) %>%
  addTiles() %>%
  addPolygons(label = ~KM2,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(KM2),
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE)
              ) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend("bottomright", pal = pal, values = ~KM2,
    title = "Municipalities extent [Km2] (DANE, 2018)",
    opacity = 1
  )
mapa
class(mun_santander$MPIO_CCDGO)
## [1] "character"
class(eva_santander$COD_MUN)
## [1] "numeric"
mun_santander$COD_MUN <-as.double(mun_santander$MPIO_CCDGO)
cacao_santander <-eva_santander%>%filter(CULTIVO=="CACAO")%>%dplyr::select(MUNICIPIO,COD_MUN,YEAR,PERIODO,TON_PROD, RENDIM)
cacao_santander
unique(cacao_santander$YEAR)
##  [1] 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018
unique(cacao_santander$PERIODO)
##  [1] "2007" "2008" "2009" "2010" "2011" "2012" "2013" "2014" "2015" "2016"
## [11] "2017" "2018"
cacao_santander%>% gather("YEAR","PERIODO","TON_PROD","RENDIM", key=variable, value=number)
summary(cacao_santander)
##   MUNICIPIO            COD_MUN           YEAR        PERIODO         
##  Length:476         Min.   :  685   Min.   :2007   Length:476        
##  Class :character   1st Qu.:68176   1st Qu.:2010   Class :character  
##  Mode  :character   Median :68322   Median :2013   Mode  :character  
##                     Mean   :60643   Mean   :2013                     
##                     3rd Qu.:68547   3rd Qu.:2016                     
##                     Max.   :68895   Max.   :2018                     
##                                                                      
##     TON_PROD          RENDIM      
##  Min.   :   0.0   Min.   :  1.00  
##  1st Qu.:  25.0   1st Qu.:  5.00  
##  Median :  77.0   Median :  6.00  
##  Mean   : 300.3   Mean   : 15.83  
##  3rd Qu.: 225.2   3rd Qu.:  8.00  
##  Max.   :6625.0   Max.   :116.00  
##                   NA's   :13
cacao_santander%>%replace(is.na(.),0)->cacao_santander2
cacao_santander %>% group_by(MUNICIPIO,COD_MUN,YEAR)%>%
  summarize(TON_PROD=sum(TON_PROD))->cacao_santander2
## `summarise()` has grouped output by 'MUNICIPIO', 'COD_MUN'. You can override using the `.groups` argument.
cacao_santander2
cacao_santander2 %>%
  group_by(COD_MUN) %>%
  gather("TON_PROD",key=variable, value=number) %>%
  unite(combi,variable,YEAR) %>%
  pivot_wider(names_from = combi, values_from = number, values_fill = 0) -> cacao_santander3
cacao_santander3
mun_santander_cacao = left_join(mun_santander,cacao_santander3, by= "COD_MUN")
mun_santander_cacao
library(RColorBrewer)
library(leaflet)
mun_santander_cacao = left_join(mun_santander, cacao_santander3, by="COD_MUN")

bins <- c(0,100,200,270,300,400,500,600,700,800,900,1000,2000)
pal <- colorBin("YlOrRd", domain = mun_santander_cacao$TON_PROD_2018, bins = bins)
 mapa <- leaflet(data = mun_santander_cacao) %>%
  addTiles() %>%
  addPolygons(label = ~TON_PROD_2018,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(TON_PROD_2018),
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE))%>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend("topright", pal = pal, values = ~TON_PROD_2017,
    title = "Producción de cacao en Santander (2018)",
    opacity = 1
  );mapa
## Warning in pal(TON_PROD_2018): Some values were outside the color scale and will
## be treated as NA
## Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette YlOrRd is 9
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette YlOrRd is 9
## Returning the palette you asked for with that many colors