Segundo Cuaderno

Introduccion

Mapa tematico para los dos mayores cultivos del departamento.

Setup

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(readr)

3. Read files related to municipios, crops and cities

list.files(pattern=c("csv"))

## [1] "boy_hort_2020.csv"                             
## [2] "boy_tuber_2020.csv"                            
## [3] "Ciudades.csv"                                  
## [4] "co.csv"                                        
## [5] "Evaluaciones_Agropecuarias_Municipales_EVA.csv"

(Tuberculos=read_csv("boy_tuber_2020.csv", show_col_types = FALSE))

## # A tibble: 121 x 4
##    Cod_Mun Municipio    Grupo                 max_prod
##      <dbl> <chr>        <chr>                    <dbl>
##  1   15740 Siachoque    Tubérculos Y Plátanos    78740
##  2   15632 Saboyá       Tubérculos Y Plátanos    63500
##  3   15861 Ventaquemada Tubérculos Y Plátanos    60144
##  4   15001 Tunja        Tubérculos Y Plátanos    59000
##  5   15764 Soracá       Tubérculos Y Plátanos    32850
##  6   15051 Arcabuco     Tubérculos Y Plátanos    23040
##  7   15646 Samacá       Tubérculos Y Plátanos    22950
##  8   15814 Toca         Tubérculos Y Plátanos    21600
##  9   15476 Motavita     Tubérculos Y Plátanos    20700
## 10   15763 Sotaquirá    Tubérculos Y Plátanos    18800
## # ... with 111 more rows

(municipios= st_read("C:/Users/USER/Desktop/Documentos Unal/Unal/Geomatica/Datos/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `C:\Users\USER\Desktop\Documentos Unal\Unal\Geomatica\Datos\15_BOYACA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS:  WGS 84

## Simple feature collection with 123 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS:  WGS 84
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                           MPIO_CRSLC
## 1          15      15001      TUNJA                                 1541
## 2          15      15022    ALMEIDA                                 1908
## 3          15      15047  AQUITANIA                                 1789
## 4          15      15051   ARCABUCO                                 1856
## 5          15      15087      BELÉN                                 1756
## 6          15      15090     BERBEO      Ordenanza 28 de Abril 9 de 1913
## 7          15      15092  BETÉITIVA                                 1754
## 8          15      15097    BOAVITA                                 1613
## 9          15      15104     BOYACÁ                                 1537
## 10         15      15106    BRICEÑO Ordenanza 14 del 25 de Julio de 1890
##    MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng  Shape_Area
## 1   119.68957      2017     BOYACÁ  0.5723744 0.009766301
## 2    57.67212      2017     BOYACÁ  0.3484692 0.004701759
## 3   942.14660      2017     BOYACÁ  1.8003115 0.076843504
## 4   137.89859      2017     BOYACÁ  0.7527090 0.011256738
## 5   163.08822      2017     BOYACÁ  0.6293489 0.013314920
## 6    58.01301      2017     BOYACÁ  0.4291743 0.004730850
## 7   101.89955      2017     BOYACÁ  0.4738184 0.008317810
## 8   145.30529      2017     BOYACÁ  0.6597822 0.011867743
## 9    48.02287      2017     BOYACÁ  0.3256140 0.003918022
## 10   64.59970      2017     BOYACÁ  0.4849753 0.005273255
##                          geometry
## 1  POLYGON ((-73.34014 5.58308...
## 2  POLYGON ((-73.36793 5.01349...
## 3  POLYGON ((-72.76242 5.63856...
## 4  POLYGON ((-73.50487 5.84347...
## 5  POLYGON ((-72.91692 6.08612...
## 6  POLYGON ((-73.0677 5.27048,...
## 7  POLYGON ((-72.81796 5.97422...
## 8  POLYGON ((-72.64907 6.43640...
## 9  POLYGON ((-73.34806 5.47411...
## 10 POLYGON ((-73.89118 5.73749...

(ciudades = read_csv("co.csv"))

## Rows: 1102 Columns: 9
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (5): city, country, iso2, admin_name, capital
## dbl (4): lat, lng, population, population_proper
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

## # A tibble: 1,102 x 9
##    city           lat   lng country  iso2  admin_name         capital population
##    <chr>        <dbl> <dbl> <chr>    <chr> <chr>              <chr>        <dbl>
##  1 Bogotá        4.61 -74.1 Colombia CO    Bogotá             primary    9464000
##  2 Medellín      6.24 -75.6 Colombia CO    Antioquia          admin      2529403
##  3 Cali          3.44 -76.5 Colombia CO    Valle del Cauca    admin      2471474
##  4 Barranquilla 11.0  -74.8 Colombia CO    Atlántico          admin      1274250
##  5 Cartagena    10.4  -75.5 Colombia CO    Bolívar            admin      1036412
##  6 Soacha        4.58 -74.2 Colombia CO    Cundinamarca       minor       995268
##  7 Palermo       2.89 -75.4 Colombia CO    Huila              minor       800000
##  8 Cúcuta        7.91 -72.5 Colombia CO    Norte de Santander admin       750000
##  9 Soledad      10.9  -74.8 Colombia CO    Atlántico          minor       698852
## 10 Pereira       4.81 -75.7 Colombia CO    Risaralda          admin       590554
## # ... with 1,092 more rows, and 1 more variable: population_proper <dbl>

sf.cities <-  st_as_sf(x = ciudades, 
                        coords = c("lng", "lat"))

sf.cities

## Simple feature collection with 1102 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -81.7006 ymin: -4.215 xmax: -67.4858 ymax: 13.3817
## CRS:           NA
## # A tibble: 1,102 x 8
##    city         country  iso2  admin_name    capital population population_prop~
##    <chr>        <chr>    <chr> <chr>         <chr>        <dbl>            <dbl>
##  1 Bogotá       Colombia CO    Bogotá        primary    9464000          7963000
##  2 Medellín     Colombia CO    Antioquia     admin      2529403          2529403
##  3 Cali         Colombia CO    Valle del Ca~ admin      2471474          2471474
##  4 Barranquilla Colombia CO    Atlántico     admin      1274250          1274250
##  5 Cartagena    Colombia CO    Bolívar       admin      1036412          1036412
##  6 Soacha       Colombia CO    Cundinamarca  minor       995268           995268
##  7 Palermo      Colombia CO    Huila         minor       800000           800000
##  8 Cúcuta       Colombia CO    Norte de San~ admin       750000           750000
##  9 Soledad      Colombia CO    Atlántico     minor       698852           342556
## 10 Pereira      Colombia CO    Risaralda     admin       590554           590554
## # ... with 1,092 more rows, and 1 more variable: geometry <POINT>

  st_crs(sf.cities) <- 4326

4. Subset data relevant to our department

sf.cities.joined <- st_join(sf.cities, municipios, join = st_within)

sf.cities.joined

## Simple feature collection with 1102 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -81.7006 ymin: -4.215 xmax: -67.4858 ymax: 13.3817
## Geodetic CRS:  WGS 84
## # A tibble: 1,102 x 17
##    city         country  iso2  admin_name    capital population population_prop~
##  * <chr>        <chr>    <chr> <chr>         <chr>        <dbl>            <dbl>
##  1 Bogotá       Colombia CO    Bogotá        primary    9464000          7963000
##  2 Medellín     Colombia CO    Antioquia     admin      2529403          2529403
##  3 Cali         Colombia CO    Valle del Ca~ admin      2471474          2471474
##  4 Barranquilla Colombia CO    Atlántico     admin      1274250          1274250
##  5 Cartagena    Colombia CO    Bolívar       admin      1036412          1036412
##  6 Soacha       Colombia CO    Cundinamarca  minor       995268           995268
##  7 Palermo      Colombia CO    Huila         minor       800000           800000
##  8 Cúcuta       Colombia CO    Norte de San~ admin       750000           750000
##  9 Soledad      Colombia CO    Atlántico     minor       698852           342556
## 10 Pereira      Colombia CO    Risaralda     admin       590554           590554
## # ... with 1,092 more rows, and 10 more variables: geometry <POINT [°]>,
## #   DPTO_CCDGO <chr>, MPIO_CCDGO <chr>, MPIO_CNMBR <chr>, MPIO_CRSLC <chr>,
## #   MPIO_NAREA <dbl>, MPIO_NANO <int>, DPTO_CNMBR <chr>, Shape_Leng <dbl>,
## #   Shape_Area <dbl>

boy.cities = dplyr::filter(sf.cities.joined, admin_name=='Boyacá')
boy.cities

## Simple feature collection with 123 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.5875 ymin: 4.8202 xmax: -72.1081 ymax: 7.0008
## Geodetic CRS:  WGS 84
## # A tibble: 123 x 17
##    city           country  iso2  admin_name capital population population_proper
##  * <chr>          <chr>    <chr> <chr>      <chr>        <dbl>             <dbl>
##  1 Tunja          Colombia CO    Boyacá     admin       209263            209263
##  2 Sogamoso       Colombia CO    Boyacá     minor       131105            131105
##  3 Duitama        Colombia CO    Boyacá     minor       126670            126670
##  4 Chiquinquirá   Colombia CO    Boyacá     minor        85274             85274
##  5 Puerto Boyacá  Colombia CO    Boyacá     minor        56930             56930
##  6 Paipa          Colombia CO    Boyacá     minor        31868             31868
##  7 Moniquirá      Colombia CO    Boyacá     minor        21402             21402
##  8 Samacá         Colombia CO    Boyacá     minor        20767             20767
##  9 Villa de Leyva Colombia CO    Boyacá     minor        18598             18598
## 10 Aquitania      Colombia CO    Boyacá     minor        16087             16087
## # ... with 113 more rows, and 10 more variables: geometry <POINT [°]>,
## #   DPTO_CCDGO <chr>, MPIO_CCDGO <chr>, MPIO_CNMBR <chr>, MPIO_CRSLC <chr>,
## #   MPIO_NAREA <dbl>, MPIO_NANO <int>, DPTO_CNMBR <chr>, Shape_Leng <dbl>,
## #   Shape_Area <dbl>

5. Make a map for the most important group of crops

library(tmap)
library(ggplot2)
library(ggrepel)
library(classInt)

class(Tuberculos$Cod_Mun)

## [1] "numeric"

class(municipios$MPIO_CCDGO)

## [1] "character"

Tuberculos$Cod_Mun = as.character(Tuberculos$Cod_Mun)

munic_tuberculos = left_join(municipios, Tuberculos, by = c("MPIO_CCDGO" = "Cod_Mun"))
munic_tuberculos

## Simple feature collection with 123 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS:  WGS 84
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                           MPIO_CRSLC
## 1          15      15001      TUNJA                                 1541
## 2          15      15022    ALMEIDA                                 1908
## 3          15      15047  AQUITANIA                                 1789
## 4          15      15051   ARCABUCO                                 1856
## 5          15      15087      BELÉN                                 1756
## 6          15      15090     BERBEO      Ordenanza 28 de Abril 9 de 1913
## 7          15      15092  BETÉITIVA                                 1754
## 8          15      15097    BOAVITA                                 1613
## 9          15      15104     BOYACÁ                                 1537
## 10         15      15106    BRICEÑO Ordenanza 14 del 25 de Julio de 1890
##    MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng  Shape_Area Municipio
## 1   119.68957      2017     BOYACÁ  0.5723744 0.009766301     Tunja
## 2    57.67212      2017     BOYACÁ  0.3484692 0.004701759   Almeida
## 3   942.14660      2017     BOYACÁ  1.8003115 0.076843504 Aquitania
## 4   137.89859      2017     BOYACÁ  0.7527090 0.011256738  Arcabuco
## 5   163.08822      2017     BOYACÁ  0.6293489 0.013314920     Belén
## 6    58.01301      2017     BOYACÁ  0.4291743 0.004730850    Berbeo
## 7   101.89955      2017     BOYACÁ  0.4738184 0.008317810 Betéitiva
## 8   145.30529      2017     BOYACÁ  0.6597822 0.011867743   Boavita
## 9    48.02287      2017     BOYACÁ  0.3256140 0.003918022    Boyacá
## 10   64.59970      2017     BOYACÁ  0.4849753 0.005273255   Briceño
##                    Grupo max_prod                       geometry
## 1  Tubérculos Y Plátanos    59000 POLYGON ((-73.34014 5.58308...
## 2  Tubérculos Y Plátanos       50 POLYGON ((-73.36793 5.01349...
## 3  Tubérculos Y Plátanos     1880 POLYGON ((-72.76242 5.63856...
## 4  Tubérculos Y Plátanos    23040 POLYGON ((-73.50487 5.84347...
## 5  Tubérculos Y Plátanos     6224 POLYGON ((-72.91692 6.08612...
## 6  Tubérculos Y Plátanos       72 POLYGON ((-73.0677 5.27048,...
## 7  Tubérculos Y Plátanos      785 POLYGON ((-72.81796 5.97422...
## 8  Tubérculos Y Plátanos      150 POLYGON ((-72.64907 6.43640...
## 9  Tubérculos Y Plátanos     4800 POLYGON ((-73.34806 5.47411...
## 10 Tubérculos Y Plátanos       32 POLYGON ((-73.89118 5.73749...

breaks <- classIntervals(munic_tuberculos$max_prod, n = 6, style = 'fisher')

## Warning in classIntervals(munic_tuberculos$max_prod, n = 6, style = "fisher"):
## var has missing values, omitted in finding classes

lab_vec <- vector(length = length(breaks$brks)-1)
rounded_breaks <- round(breaks$brks,2)
lab_vec[1] <- paste0('[', rounded_breaks[1],' - ', rounded_breaks[2],']')
for(i in 2:(length(breaks$brks) - 1)){
  lab_vec[i] <- paste0('(',rounded_breaks[i], ' - ', rounded_breaks[i+1], ']')
}

munic_tuberculos <-  munic_tuberculos %>%
  mutate(faktor_class = factor(cut(max_prod, breaks$brks, include.lowest = T), labels = lab_vec))
munic_tuberculos$Produccion = munic_tuberculos$faktor_class

munic_tuberculos$mid <- sf::st_centroid(munic_tuberculos$geometry)
LONG = st_coordinates(munic_tuberculos$mid)[,1]
LAT = st_coordinates(munic_tuberculos$mid)[,2]
ggplot(data = munic_tuberculos) +
   geom_sf(aes(fill = Produccion)) +
   geom_label_repel(aes(x = LONG, y = LAT, label = MPIO_CNMBR), 
                    label.padding =     unit(0.05,"lines"),  
                    label.r = unit(0.025, "lines"),
                    label.size = 0.05)

## Warning: ggrepel: 111 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps