Brazil Municipalities SPATIAL DATA ANALYSIS
BLESSING NEHOHWA
2022-10-21
##install.packages('tidyverse')
##install.packages('sf')
##install.packages('scales')
##install.packages('ggplot2')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.
1a). Open municipalities and filter for Sp state and print selected
data within the
console
municipalities<- read_sf("C:/Users/beemu/OneDrive/SDA/municipalities/municipalities.shp" )
municipalities_SP<-municipalities%>% filter(UF=="SP")
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 × 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…1b). Generating municipalities centroids and plotting them
municipalities_centroids<-municipalities_SP%>%st_centroid()
municipalities_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.
1c).Dissolve SP State polygons to create SP state Boarder
SP_boarder<-municipalities_SP%>%st_union()
SP_boarder%>% ggplot()+geom_sf()
1d). Plotting SP centroids and SP boarder together
SP_boarder%>%ggplot() +geom_sf() +
geom_sf(data=municipalities_centroids) +
theme_classic()
Question 2
What is the mean Human Development Index of Municipalities in each state of Brazil? HDI is the variable from Municipalities data set.
2a).Reading The table
br_states <-municipalities[,]
print(br_states)## Simple feature collection with 5840 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -73.95055 ymin: -33.74511 xmax: -32.38177 ymax: 5.272156
## Geodetic CRS: WGS 84
## # A tibble: 5,840 × 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, -4…
## 2 1507953 TAILANDIA PA 85468 0.588 5405 ((-48.92547 -3.406675, -…
## 3 1500206 ACARA PA 53787 0.506 4389 ((-48.49495 -2.558588, -…
## 4 1508001 TOME ACU PA 57914 0.586 4765 ((-48.53014 -3.195302, -…
## 5 1501808 BREVES PA 94779 0.503 3608 ((-49.9747 -1.30028, -50…
## 6 1502806 CURRALINHO PA 29838 0.502 2270 ((-50.25328 -1.805703, -…
## 7 1507706 SAO SEBASTIA… PA 23696 0.558 3116 ((-49.58361 -1.730366, -…
## 8 1507706 SAO SEBASTIA… PA 23696 0.558 3116 ((-49.58429 -1.76721, -4…
## 9 1505205 OEIRAS DO PA… PA 29402 0.507 3933 ((-49.72789 -1.846166, -…
## 10 1504000 LIMOEIRO DO … PA 25846 0.541 3767 ((-49.70602 -2.068686, -…
## # … with 5,830 more rows2b). Calculating the mean HDI by state
br_states <-municipalities %>% group_by(UF) %>% summarize(IDHM_mean=mean(IDHM_10))
br_states2c). Plotting the result (Mean)
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”
3a). Import indigenous park shapefile and plot it. The projections are the same
indigenous<-read_sf("C:/Users/beemu/OneDrive/SDA/BC250_Terra_Indigena_A/BC250_Terra_Indigena_A.shp") %>% st_transform(4326)
br_states %>% ggplot() +
geom_sf() +
geom_sf(data=indigenous,fill="red")
3b). Creating 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)3c). Plot the shapefiles to see the overlap
Gaucha%>% ggplot() +
geom_sf(fill="red") +
geom_sf(data=Xingu,fill="blue",alpha=0.5)
3d).Creating shapefile for intersection and plotting it Attribute variables are assumed to be spatially constant throughout all geometries changed fill colour to Orange
intersection <- Gaucha %>% st_intersection(Xingu)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesintersection %>% ggplot() + geom_sf(fill="black")
3e). 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")
3e.) The area of the intersection
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.
4a). Import Housing Shapefile
Housing <- read_sf("C:/Users/beemu/OneDrive/SDA/MCMV_new/MCMV_new.shp")4b). Select housing in AC state and plot the result.
Housing_AC <- Housing %>% filter(UF=="AC")
br_states %>% filter(UF=="AC") %>% ggplot() +
geom_sf(fill="light blue") + geom_sf(data = Housing_AC)
4c). Calculating 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]4d). Calculating buffer around housing points
Housing_AC %>% st_transform(29189) %>%
st_buffer(20000) %>%
ggplot() + geom_sf(fill="dark green")
4e). Plotting it all Together
br_states %>% filter(UF=="AC") %>% ggplot() +
geom_sf(fill="light blue") +
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.
5a). 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 × 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>5b). Calculating 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.5c). Select municipality with most Housing units
mun_Housing_units %>% arrange(-UH) %>% slice(1) %>% pull(NOME)## [1] "RIO DE JANEIRO"5d) Plotting 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 blue",trans="log") +
theme_classic()## Warning: Transformation introduced infinite values in discrete y-axis
NB: Transformation introduced infinite values in discrete y-axis
##EXTRA: CREATE A VORONOI POLYGON FOR SP STATE
Extra 1a). Installing and Loading needed libraries
#install.packages("dismo")
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)## Loading required package: raster## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:ggsn':
##
## scalebar## The following object is masked from 'package:dplyr':
##
## selectExtra 1b). Creating and viewing 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(fill="white")
Extra 1c). Clipping Voronoi to SP State
sp_voronoi_clip <- sp_voronoi %>% st_intersection(st_union(SP_boarder))## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriessp_voronoi_clip %>% ggplot() + geom_sf()
Warning: attribute variables are assumed to be spatially constant
throughout all geometries.
EXtra 2: Count points inside a polygon
Perform operation and view resultsmun_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 = "purple") +
theme_classic()