Question 1

Produce a map showing the centroids of each municipality in just the state of Säo Paulo, and add the outer boundary of Säo Paulo state.

  1. Open municipalities and filter for SP state
municipalites <- read_sf("/Users/ogeohia/Downloads/municipalities (1) 2/municipalities.shp")

municipalites_SP <- municipalites %>% filter(UF=="SP")
  1. Generating municipalities centroids and plotting them
municipalites_centroids <- municipalites_SP %>% st_centroid()

municipalites_centroids %>% ggplot() + geom_sf()

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

  1. Dissolve SP state polygons to create SP state border
SP_border <- municipalites_SP %>% st_union()

SP_border %>% ggplot() + geom_sf()

  1. Plot SP centroids and SP border together
SP_border %>% ggplot() + geom_sf() + 
  geom_sf(data=municipalites_centroids) + 
  theme_classic()

Question 2

What is the mean Human Development Index of municipalities in each state of Brazil?

  1. View municipalities data table
head(municipalites)
## Simple feature collection with 6 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -51.21642 ymin: -3.414883 xmax: -47.93154 ymax: -0.560798
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 7
##   COD_MUN NOME       UF    POP_201 IDHM_10 PIB_PER                      geometry
##     <int> <chr>      <chr>   <dbl>   <dbl>   <dbl>                 <POLYGON [°]>
## 1 1504703 MOJU       PA      72597   0.547    3894 ((-48.8278 -2.601164, -48.83…
## 2 1507953 TAILANDIA  PA      85468   0.588    5405 ((-48.92547 -3.406675, -48.9…
## 3 1500206 ACARA      PA      53787   0.506    4389 ((-48.49495 -2.558588, -48.5…
## 4 1508001 TOME ACU   PA      57914   0.586    4765 ((-48.53014 -3.195302, -48.4…
## 5 1501808 BREVES     PA      94779   0.503    3608 ((-49.9747 -1.30028, -50.017…
## 6 1502806 CURRALINHO PA      29838   0.502    2270 ((-50.25328 -1.805703, -50.3…
  1. Calculate mean HDI by state and view result data table
br_states <- municipalites %>% group_by(UF) %>% summarize(IDHM_mean=mean(IDHM_10))
head(br_states)
## Simple feature collection with 6 features and 2 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -73.95055 ymin: -18.3494 xmax: -35.15229 ymax: 4.43706
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 3
##   UF    IDHM_mean                                                       geometry
##   <chr>     <dbl>                                                 <GEOMETRY [°]>
## 1 AC        0.586 POLYGON ((-69.34267 -10.95212, -69.19872 -10.94806, -69.01297…
## 2 AL        0.564 POLYGON ((-35.90323 -9.869403, -35.85563 -9.789107, -35.79578…
## 3 AM        0.565 POLYGON ((-60.75412 -0.841936, -60.76055 -0.761119, -60.975 -…
## 4 AP        0.658 MULTIPOLYGON (((-50.5029 2.099758, -50.52607 2.009225, -50.48…
## 5 BA        0.596 MULTIPOLYGON (((-46.04805 -13.30876, -46.1336 -13.37479, -46.…
## 6 CE        0.617 POLYGON ((-38.40454 -6.056894, -38.33925 -6.077296, -38.27457…
  1. Plotting the result
br_states %>% ggplot() + geom_sf(aes(fill=IDHM_mean)) +
  theme_classic() + scale_fill_continuous(low="red", high="green")

Question 3

Produce a polygon/shapefile mapping the area of the municipality ‘Gaucha do Norte’ that is in the indigenous territory “Parque do Xingu”.

  1. Importing and plotting the indigenous park shapefile
indigenous <- read_sf("/Users/ogeohia/Downloads/BC250_Terra_Indigena_A/BC250_Terra_Indigena_A.shp") %>% st_transform(4326)

br_states %>% ggplot() + geom_sf() + geom_sf(data=indigenous, fill="red")

  1. Creating a shapefile for Gaucha do Norte in Parque do Xingu
Xingu <- indigenous %>% filter(nome=="Parque do Xingu") %>% st_transform(4326)
Gaucha <- municipalites %>% filter(NOME=="GAUCHA DO NORTE") %>% st_transform(4326)
  1. Plotting the shapefiles to see the overlap
Gaucha %>% ggplot() + geom_sf(fill="red") + geom_sf(data=Xingu, fill="blue", alpha=0.5)

  1. Creating shapefile for intersection and plotting it
intersection <- Gaucha %>% st_intersection(Xingu)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
intersection %>% ggplot() + geom_sf(fill="dark green")

  1. Plotting all together
Gaucha %>% ggplot() + geom_sf(fill="red", alpha=0.5) + geom_sf(data=Xingu, fill="blue", alpha=0.5) + geom_sf(data=intersection, fill="black")

You can calculate the intersection area

st_area(intersection)
## 8109249396 [m^2]

Question 4

In the state of Acre (AC), which two social housing (MCMV) projects are closest to each other? Create a 10km buffer around each housing project.

  1. Import housing shapefile
Housing <- read_sf("/Users/ogeohia/Downloads/MCMV_new/MCMV_new.shp")
  1. Selecting housing in AC state and plotting the result
Housing_AC <- Housing %>% filter(UF=="AC")

br_states %>% filter(UF == "AC") %>% ggplot() + geom_sf() + geom_sf(data = Housing_AC)

  1. Calculating the distance between housing pints and view results
distance <- Housing_AC %>% 
  st_transform(29189) %>% 
  st_distance() %>% 
  as.data.frame()
head(distance)
##               1             2             3             4            5
## 1      0.00 [m]  24804.69 [m] 242315.19 [m] 287859.10 [m] 422339.2 [m]
## 2  24804.69 [m]      0.00 [m] 218072.03 [m] 263435.24 [m] 398549.5 [m]
## 3 242315.19 [m] 218072.03 [m]      0.00 [m]  46436.01 [m] 181680.9 [m]
## 4 287859.10 [m] 263435.24 [m]  46436.01 [m]      0.00 [m] 139918.1 [m]
## 5 422339.25 [m] 398549.47 [m] 181680.92 [m] 139918.06 [m]      0.0 [m]
## 6 592605.19 [m] 572944.86 [m] 387497.33 [m] 361461.06 [m] 248657.5 [m]
##              6             7            8            9           10
## 1 592605.2 [m] 492622.07 [m] 596084.1 [m] 617807.7 [m] 704344.7 [m]
## 2 572944.9 [m] 468877.05 [m] 573086.5 [m] 594959.9 [m] 680777.6 [m]
## 3 387497.3 [m] 251979.33 [m] 359890.3 [m] 382594.3 [m] 464241.1 [m]
## 4 361461.1 [m] 209540.53 [m] 320046.9 [m] 343176.1 [m] 421487.9 [m]
## 5 248657.5 [m]  70363.19 [m] 180507.2 [m] 203859.1 [m] 282581.5 [m]
## 6      0.0 [m] 216294.81 [m] 158054.3 [m] 153476.8 [m] 240517.0 [m]
  1. Generating buffer around housing points
Housing_AC %>% st_transform(29189) %>% st_buffer(20000) %>% ggplot() + geom_sf(fill="dark green")

  1. Plotting it all together
br_states %>% filter(UF == "AC") %>% ggplot() + 
  geom_sf() + 
  geom_sf(data = (Housing_AC %>% st_transform(29189) %>% st_buffer(20000)), fill = "dark green") +
  geom_sf(data = Housing_AC, color = "red")

Question 5

Across Brazil, which municipalities have the lowest and highest number of MCMV housing units (UH) in its territory? Create a map of the distribution of total housing units by municipality.

  1. Performing spatial join between municipalities and housing
mun_Housing_units <- municipalites %>% st_join(Housing)
  1. Calculating total of Housing units in each municipality
mun_Housing_units <- mun_Housing_units %>% 
  group_by(COD_MUN, NOME) %>% 
  summarize(UH=sum(UH,na.rm=T)) %>% 
  ungroup()
## `summarise()` has grouped output by 'COD_MUN'. You can override using the
## `.groups` argument.
  1. Selecting municipality with the most Housing units
mun_Housing_units %>% arrange(-UH) %>% slice(1) %>% pull(NOME)
## [1] "RIO DE JANEIRO"
  1. Plotting the distribution of total housing units by municipality
mun_Housing_units %>% ggplot() +
  geom_sf(aes(fill=UH), col=NA) +
  scale_fill_gradient2(low="#ccece6", high="dark green", trans="log") +
  theme_classic()
## Warning in scale_fill_gradient2(low = "#ccece6", high = "dark green", trans =
## "log"): log-2.718282 transformation introduced infinite values.