SDA WEEK 3 FORMATIVE ASSESSMENT
Falaye Bankole ojo
8/19/2022
First task
Produce a map showing centroids of each municipality in just the state of Sao Paulo, and add the outer boundary of Sao Paulo state.
Open municipalities and filter for SP.
municipalities <- read_sf("OGIS Brazilian Munipalities/municipalities.shp")
municipalities_sp <- municipalities %>% filter(UF=="SP")Prints the selected data within the console.
print(municipalities_sp[1:10,])## Simple feature collection with 10 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -52.65639 ymin: -25.31024 xmax: -47.73479 ymax: -22.1845
## Geodetic CRS: WGS 84
## # A tibble: 10 x 7
## COD_MUN NOME UF POP_201 IDHM_10 PIB_PER geometry
## <int> <chr> <chr> <dbl> <dbl> <dbl> <POLYGON [°]>
## 1 3509908 CANANEIA SP 12216 0.72 9201 ((-48.09754 -25.31023, -~
## 2 3509908 CANANEIA SP 12216 0.72 9201 ((-47.91403 -25.16738, -~
## 3 3509908 CANANEIA SP 12216 0.72 9201 ((-47.86343 -25.12852, -~
## 4 3505351 BARRA DO CHA~ SP 5305 0.66 7305 ((-49.25004 -24.44493, -~
## 5 3522653 ITAPIRAPUA P~ SP 3926 0.661 6822 ((-49.20744 -24.70042, -~
## 6 3523206 ITARARE SP 48143 0.703 12292 ((-49.25389 -24.31763, -~
## 7 3507159 BOM SUCESSO ~ SP 3623 0.66 8507 ((-49.25192 -24.40347, -~
## 8 3554300 TEODORO SAMP~ SP 21595 0.741 12599 ((-52.28946 -22.62776, -~
## 9 3504008 ASSIS SP 96336 0.805 14272 ((-50.5712 -22.68999, -5~
## 10 3528809 MARACAI SP 13382 0.771 32164 ((-50.5138 -22.6032, -50~Generate municipalities centroids and plot them
municipalities_centroids <- municipalities_sp %>% st_centroid()
municipalities_centroids %>% ggplot() + geom_sf()
Dissolve SP state polygons to create SP state border
SP_border <- municipalities_sp %>% st_union()
SP_border %>% ggplot() + geom_sf()
Plot SP centroids and border together
SP_border %>% ggplot() + geom_sf() + geom_sf(data=municipalities_centroids) + theme_classic()
Second task
What is the mean Human Development Index of municipalities in each state of brazil?
1. Remember that HDI is a variable from a municipalities data set. Open table to see it.
print(municipalities[1:4, 3:5])## Simple feature collection with 4 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -49.62303 ymin: -3.414883 xmax: -47.93154 ymax: -1.469574
## Geodetic CRS: WGS 84
## # A tibble: 4 x 4
## UF POP_201 IDHM_10 geometry
## <chr> <dbl> <dbl> <POLYGON [°]>
## 1 PA 72597 0.547 ((-48.8278 -2.601164, -48.83234 -2.635788, -48.93922 -2~
## 2 PA 85468 0.588 ((-48.92547 -3.406675, -48.96627 -3.212913, -48.91762 -~
## 3 PA 53787 0.506 ((-48.49495 -2.558588, -48.567 -2.482825, -48.66086 -2.~
## 4 PA 57914 0.586 ((-48.53014 -3.195302, -48.49342 -3.0943, -48.5534 -2.9~2. Calculate mean HDI by state.
br_states <- municipalities %>% group_by(UF) %>% summarize(IDHM_mean=mean(IDHM_10))
print(br_states[1:5,])## Simple feature collection with 5 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: 5 x 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.~3. Plot the result
br_states %>% ggplot() +
geom_sf(aes(fill=IDHM_mean)) +
theme_classic() +
scale_fill_continuous(low="red", high="green")
Third task
Produce a polygon/shapefile mapping the area of the municipality ‘Gaucha do norte’ that is in the indigenous territory
“Parque do Xingu”
1. Import the indigenous park shapefile and plot
indigenous <- read_sf("Indegenous Groups/BC250_Terra_Indigena_A.shp") %>% st_transform(4326)
br_states %>% ggplot() +
geom_sf() +
geom_sf(data=indigenous, fill="red")
2. Create a shapefile for Gaucha do Norte e Xingu.
Xingu <- indigenous %>% filter(nome=="Parque do Xingu") %>% st_transform(4326)
Gaucha <- municipalities %>% filter(NOME=="GAUCHA DO NORTE") %>% st_transform(4326)3. Plot the shapefiles to see the overlap.
Gaucha %>% ggplot() +
geom_sf(fill="red") +
geom_sf(data=Xingu, fill="blue", alpha=0.5)
4. Create shapefile for intersection and plot it.
intersection <- Gaucha %>% st_intersection(Xingu)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesintersection %>% ggplot() + geom_sf(fill="dark green")
5. Plot 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, for example.
st_area(intersection)## 8109249396 [m^2]Fourth task
1. Import housing shapefile
Housing <- read_sf("Housing/MCMV_new.shp")2. Select housing in AC state and plot the result.
Housing_AC <- Housing %>% filter(UF=="AC")
br_states %>% filter(UF=="AC") %>% ggplot() +
geom_sf() + geom_sf(data = Housing_AC)
3. Calculate distance between housing points and view results
distance <- Housing_AC %>%
st_transform(29189) %>%
st_distance() %>%
as.data.frame()
print(distance[1:5,])## 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 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]4. Calculate buffer around housing points
Housing_AC %>% st_transform(29189) %>%
st_buffer(20000) %>%
ggplot() + geom_sf(fill="dark green")
5. Plot 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")
Fifth task
Across Brazil, which municipality has the highest number of MCMV housing units (HU) in its teritory? create the a map of the distribution of total housing units by municipality.
1. spatial join between municipalities and housing.
mun_Housing_units <-municipalities %>% st_join(Housing)
print(mun_Housing_units[1:5,])## Simple feature collection with 5 features and 11 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: 5 x 12
## COD_MUN NOME UF.x POP_201 IDHM_10 PIB_PER geometry XCOORD
## <int> <chr> <chr> <dbl> <dbl> <dbl> <POLYGON [°]> <dbl>
## 1 1504703 MOJU PA 72597 0.547 3894 ((-48.8278 -2.601164, -4~ -48.8
## 2 1507953 TAILAN~ PA 85468 0.588 5405 ((-48.92547 -3.406675, -~ -48.9
## 3 1500206 ACARA PA 53787 0.506 4389 ((-48.49495 -2.558588, -~ NA
## 4 1508001 TOME A~ PA 57914 0.586 4765 ((-48.53014 -3.195302, -~ -48.2
## 5 1501808 BREVES PA 94779 0.503 3608 ((-49.9747 -1.30028, -50~ -50.5
## # ... with 4 more variables: YCOORD <dbl>, UF.y <chr>, UH <dbl>,
## # Project_ID <int>2. Calculate total of Housing units in each municipality
mun_Housing_units <- mun_Housing_units %>%
group_by(COD_MUN,NOME) %>%
summarise(UH=sum(UH,na.rm = T)) %>%
ungroup()## `summarise()` has grouped output by 'COD_MUN'. You can override using the
## `.groups` argument.3. Select municipality with most Housing units
mun_Housing_units %>% arrange(-UH) %>% slice(1) %>% pull(NOME)## [1] "RIO DE JANEIRO"4. Plot the distribution of total Housing units by municipality
mun_Housing_units %>% ggplot() +
geom_sf(aes(fill=UH), col=NA) +
scale_fill_gradient(low="#ccece6",high="dark red",trans="log") +
theme_classic()## Warning: Transformation introduced infinite values in discrete y-axis
EXTRA: Create a vpronoi polygon for SP state.
1. Load needed libraries.
library(deldir)## deldir 1.0-6 Nickname: "Mendacious Cosmonaut"##
## The syntax of deldir() has had an important change.
## The arguments have been re-ordered (the first three
## are now "x, y, z") and some arguments have been
## eliminated. The handling of the z ("tags")
## argument has been improved.
##
## The "dummy points" facility has been removed.
## This facility was a historical artefact, was really
## of no use to anyone, and had hung around much too
## long. Since there are no longer any "dummy points",
## the structure of the value returned by deldir() has
## changed slightly. The arguments of plot.deldir()
## have been adjusted accordingly; e.g. the character
## string "wpoints" ("which points") has been
## replaced by the logical scalar "showpoints".
## The user should consult the help files.library(dismo)## Warning: package 'dismo' was built under R version 4.1.3## Loading required package: raster## Warning: package 'raster' was built under R version 4.1.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.1.3##
## Attaching package: 'raster'## The following object is masked from 'package:ggsn':
##
## scalebar## The following object is masked from 'package:dplyr':
##
## select2. Create and view voronoi polygon
sp_voronoi <- municipalities_centroids %>%
as("Spatial") %>%
voronoi() %>%
st_as_sf()## Warning in proj4string(x): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.htmlsp_voronoi %>% ggplot() +
geom_sf()
3. Clip voronoi to SP state border (clip example). Remember to check projection!
sp_voronoi_clip <- sp_voronoi %>% st_intersection(st_union(SP_border))## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriessp_voronoi_clip %>% ggplot() + geom_sf()
EXTRA 2: count points inside a polygon
1. Perform operation and view results
mun_housin <- municipalities %>% st_join(Housing) %>%
group_by(COD_MUN) %>%
count()
mun_housin %>% ggplot() +
geom_sf(aes(fill=n),col=NA) +
scale_fill_gradient(low = "#ccece6", high = "dark green") +
theme_classic()