Map making in R (GIS)
Dr Gurpreet Singh, Dr Biju Soman
18/07/2019
R can be used for making and arranging maps in GIS using packages namely tidyverse, sf, ggplot2, ggrepel ``` As a first step, clear the environment using following commands
#To see the list of files loaded in environment
ls() ## character(0)#To remove all the files from the environment
rm(list=ls())
#To confirm that the environment is clean
ls()## character(0)Set your working directory to the location where files are stored in the computer. To do so, click on three dot sign present in the panel of files interface (right lower panel in R studio) -> Go to the file location -> click on more button in the file panel and click on "set as working directory". Alternatively, working directory can be set using command:-#setwd("C:/Users/Gurpreet/Desktop/r scripts")
#To check/ confirm the location of working directory
#getwd()
#Note: The location of wd can also be seen on the panel of console.Load the libraries - tidyverse and sf. In case of first time use, first install the packages/ libraries and then load using command:-
#To install, use command install.packages(c("sf", "tidyverse", "tmap"))
#To load the libraries
library(tidyverse)## -- Attaching packages ----------------------------------------- tidyverse 1.2.1 --## v ggplot2 3.1.1 v purrr 0.3.2
## v tibble 2.1.1 v dplyr 0.8.0.1
## v tidyr 0.8.3 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0## -- Conflicts -------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library (tmap)Reading rds file. Use readRDS command to read file from the wd. Alternatively, in case the file has to be read from other location in the computer, file.choose command can be used #simp_sf <- readRDS("simp_sf.rds")
#Alternate
simp_sf = readRDS(file.choose())Its a good habit to glimpse through the file to understand the variables#to see the summary of the file
summary(simp_sf)## state_ut abb type data_year
## Length:36 AN : 1 Length:36 Length:36
## Class :character AP : 1 Class :character Class :character
## Mode :character AR : 1 Mode :character Mode :character
## AS : 1
## BR : 1
## CH : 1
## (Other):30
## comparable_economy nominal_gdp_usd nominal_gdp_inr
## Length:36 Min. :6.000e+07 Min. :4.070e+09
## Class :character 1st Qu.:4.100e+09 1st Qu.:2.705e+11
## Mode :character Median :4.650e+10 Median :3.035e+12
## Mean :8.051e+10 Mean :5.207e+12
## 3rd Qu.:1.268e+11 3rd Qu.:8.295e+12
## Max. :4.300e+11 Max. :2.796e+13
##
## pop_2011 decadal_growth rural_pop
## Min. : 64429 Min. :-0.0050 Min. : 14121
## 1st Qu.: 1438945 1st Qu.: 0.1385 1st Qu.: 545820
## Median : 21070511 Median : 0.1829 Median : 12833156
## Mean : 33624330 Mean : 0.1900 Mean : 23142433
## 3rd Qu.: 52136006 3rd Qu.: 0.2238 3rd Qu.: 34820100
## Max. :199812341 Max. : 0.5550 Max. :155111022
##
## urban_pop area_km2 density_km2 sex_ratio
## Min. : 50308 Min. : 32 Min. : 17.0 Min. : 618.0
## 1st Qu.: 665287 1st Qu.: 9927 1st Qu.: 174.8 1st Qu.: 904.2
## Median : 5162647 Median : 54578 Median : 334.5 Median : 946.5
## Mean :10471420 Mean : 91389 Mean : 1069.9 Mean : 931.6
## 3rd Qu.:16032607 3rd Qu.:140320 3rd Qu.: 730.5 3rd Qu.: 975.8
## Max. :50827531 Max. :342239 Max. :11297.0 Max. :1084.0
##
## region per_capita_gdp_inr per_capita_gdp_usd
## Length:36 Min. : 43596 Min. : 658.7
## Class :character 1st Qu.: 97191 1st Qu.:1466.8
## Mode :character Median :150970 Median :2194.5
## Mean :175026 Mean :2706.5
## 3rd Qu.:241112 3rd Qu.:3715.3
## Max. :482945 Max. :7546.0
##
## geometry
## MULTIPOLYGON :36
## epsg:4326 : 0
## +proj=long...: 0
##
##
##
## #to see the top 6 rows in the file (by default setting)
head(simp_sf)## Simple feature collection with 6 features and 17 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 72.57426 ymin: 12.62373 xmax: 84.71323 ymax: 37.04546
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## state_ut abb type data_year comparable_economy nominal_gdp_usd
## 1 Andhra Pradesh AP State 2018–19 Hungary 1.36e+11
## 2 Jammu & Kashmir JK State 2018–19 Iceland 2.50e+10
## 3 Himachal Pradesh HP State 2018–19 Iceland 2.50e+10
## 4 Punjab PB State 2018–19 Sri Lanka 8.00e+10
## 5 Uttarakhand UT State 2018–19 Turkmenistan 4.00e+10
## 6 Haryana HR State 2018–19 Slovakia 1.10e+11
## nominal_gdp_inr pop_2011 decadal_growth rural_pop urban_pop area_km2
## 1 8.70e+12 49386799 0.111 34776389 14610410 162968
## 2 1.57e+12 12548926 0.237 9134820 3414106 222236
## 3 1.52e+12 6864602 0.128 6167805 688704 55673
## 4 5.18e+12 27704236 0.137 17316800 10387436 50362
## 5 2.58e+12 10116752 0.192 7025583 3091169 53483
## 6 6.87e+12 25353081 0.199 16531493 8821588 44212
## density_km2 sex_ratio region per_capita_gdp_inr per_capita_gdp_usd
## 1 303 993 Southern 176160.4 2753.772
## 2 57 883 Northern 125110.3 1992.202
## 3 123 974 Northern 221425.8 3641.872
## 4 550 893 Northern 186975.0 2887.645
## 5 189 963 Northern 255022.6 3953.838
## 6 573 877 Northern 270973.0 4338.723
## geometry
## 1 MULTIPOLYGON (((80.26915 13...
## 2 MULTIPOLYGON (((78.4342 32....
## 3 MULTIPOLYGON (((78.4342 32....
## 4 MULTIPOLYGON (((75.35152 32...
## 5 MULTIPOLYGON (((78.9706 31....
## 6 MULTIPOLYGON (((76.83514 30...#to see only top 3 rows
head(simp_sf,3) ## Simple feature collection with 3 features and 17 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 72.57426 ymin: 12.62373 xmax: 84.71323 ymax: 37.04546
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## state_ut abb type data_year comparable_economy nominal_gdp_usd
## 1 Andhra Pradesh AP State 2018–19 Hungary 1.36e+11
## 2 Jammu & Kashmir JK State 2018–19 Iceland 2.50e+10
## 3 Himachal Pradesh HP State 2018–19 Iceland 2.50e+10
## nominal_gdp_inr pop_2011 decadal_growth rural_pop urban_pop area_km2
## 1 8.70e+12 49386799 0.111 34776389 14610410 162968
## 2 1.57e+12 12548926 0.237 9134820 3414106 222236
## 3 1.52e+12 6864602 0.128 6167805 688704 55673
## density_km2 sex_ratio region per_capita_gdp_inr per_capita_gdp_usd
## 1 303 993 Southern 176160.4 2753.772
## 2 57 883 Northern 125110.3 1992.202
## 3 123 974 Northern 221425.8 3641.872
## geometry
## 1 MULTIPOLYGON (((80.26915 13...
## 2 MULTIPOLYGON (((78.4342 32....
## 3 MULTIPOLYGON (((78.4342 32....#to see names of the variables
names(simp_sf)## [1] "state_ut" "abb" "type"
## [4] "data_year" "comparable_economy" "nominal_gdp_usd"
## [7] "nominal_gdp_inr" "pop_2011" "decadal_growth"
## [10] "rural_pop" "urban_pop" "area_km2"
## [13] "density_km2" "sex_ratio" "region"
## [16] "per_capita_gdp_inr" "per_capita_gdp_usd" "geometry"#to see dimensions of the file
dim(simp_sf)## [1] 36 18Map making. # to plot a selected variable/ column
plot(simp_sf['pop_2011'])
#to plot only geometry of the map
plot(simp_sf[0])
#alternative method
plot(st_geometry(simp_sf))
Tip...shortcut to make a pipe command '%>%' is "cntrl+shift+M"
Creation of maps with tmap ```
library(tmap)
head(simp_sf)## Simple feature collection with 6 features and 17 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 72.57426 ymin: 12.62373 xmax: 84.71323 ymax: 37.04546
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## state_ut abb type data_year comparable_economy nominal_gdp_usd
## 1 Andhra Pradesh AP State 2018–19 Hungary 1.36e+11
## 2 Jammu & Kashmir JK State 2018–19 Iceland 2.50e+10
## 3 Himachal Pradesh HP State 2018–19 Iceland 2.50e+10
## 4 Punjab PB State 2018–19 Sri Lanka 8.00e+10
## 5 Uttarakhand UT State 2018–19 Turkmenistan 4.00e+10
## 6 Haryana HR State 2018–19 Slovakia 1.10e+11
## nominal_gdp_inr pop_2011 decadal_growth rural_pop urban_pop area_km2
## 1 8.70e+12 49386799 0.111 34776389 14610410 162968
## 2 1.57e+12 12548926 0.237 9134820 3414106 222236
## 3 1.52e+12 6864602 0.128 6167805 688704 55673
## 4 5.18e+12 27704236 0.137 17316800 10387436 50362
## 5 2.58e+12 10116752 0.192 7025583 3091169 53483
## 6 6.87e+12 25353081 0.199 16531493 8821588 44212
## density_km2 sex_ratio region per_capita_gdp_inr per_capita_gdp_usd
## 1 303 993 Southern 176160.4 2753.772
## 2 57 883 Northern 125110.3 1992.202
## 3 123 974 Northern 221425.8 3641.872
## 4 550 893 Northern 186975.0 2887.645
## 5 189 963 Northern 255022.6 3953.838
## 6 573 877 Northern 270973.0 4338.723
## geometry
## 1 MULTIPOLYGON (((80.26915 13...
## 2 MULTIPOLYGON (((78.4342 32....
## 3 MULTIPOLYGON (((78.4342 32....
## 4 MULTIPOLYGON (((75.35152 32...
## 5 MULTIPOLYGON (((78.9706 31....
## 6 MULTIPOLYGON (((76.83514 30...#filter command can be used to use subset of a map without creating separate file
#%in% command read as "includes"
#tm_fill command to tell how to fill and what is the title of the legend
#tm_borders tells about line width "lwd"
#tm_text command used to abbreviate names and set size of the label
#tm_layout gives main title
#tm_credits to tell the source details
#to go to next line "\n" command is used
simp_sf %>%
filter(!state_ut %in% c("Andaman & Nicobar Islands", "Lakshadweep")) %>%
tm_shape() +
tm_fill(col = "pop_2011", title = "No. People") +
tm_borders(lwd = 0.5) +
tm_text("abb", size = 0.5) +
tm_style("gray") +
tm_layout(
main.title = "Population (2011)",
main.title.position = c("center"),
main.title.size = 1,
legend.position = c("right", "bottom")
) +
tm_credits("Data:\n2011 Census", position = c("left", "bottom"))
# Creation and arrangement of maps
# functions of tmap_arrange
head(simp_sf)## Simple feature collection with 6 features and 17 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 72.57426 ymin: 12.62373 xmax: 84.71323 ymax: 37.04546
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## state_ut abb type data_year comparable_economy nominal_gdp_usd
## 1 Andhra Pradesh AP State 2018–19 Hungary 1.36e+11
## 2 Jammu & Kashmir JK State 2018–19 Iceland 2.50e+10
## 3 Himachal Pradesh HP State 2018–19 Iceland 2.50e+10
## 4 Punjab PB State 2018–19 Sri Lanka 8.00e+10
## 5 Uttarakhand UT State 2018–19 Turkmenistan 4.00e+10
## 6 Haryana HR State 2018–19 Slovakia 1.10e+11
## nominal_gdp_inr pop_2011 decadal_growth rural_pop urban_pop area_km2
## 1 8.70e+12 49386799 0.111 34776389 14610410 162968
## 2 1.57e+12 12548926 0.237 9134820 3414106 222236
## 3 1.52e+12 6864602 0.128 6167805 688704 55673
## 4 5.18e+12 27704236 0.137 17316800 10387436 50362
## 5 2.58e+12 10116752 0.192 7025583 3091169 53483
## 6 6.87e+12 25353081 0.199 16531493 8821588 44212
## density_km2 sex_ratio region per_capita_gdp_inr per_capita_gdp_usd
## 1 303 993 Southern 176160.4 2753.772
## 2 57 883 Northern 125110.3 1992.202
## 3 123 974 Northern 221425.8 3641.872
## 4 550 893 Northern 186975.0 2887.645
## 5 189 963 Northern 255022.6 3953.838
## 6 573 877 Northern 270973.0 4338.723
## geometry
## 1 MULTIPOLYGON (((80.26915 13...
## 2 MULTIPOLYGON (((78.4342 32....
## 3 MULTIPOLYGON (((78.4342 32....
## 4 MULTIPOLYGON (((75.35152 32...
## 5 MULTIPOLYGON (((78.9706 31....
## 6 MULTIPOLYGON (((76.83514 30... states_sf <- simp_sf %>% filter(!type == "Union Territory")#creating a growth map
growth <- tm_shape(states_sf) +
tm_fill(col = "decadal_growth", title = "Percentage") +
tm_borders(lwd = 0.5) +
tm_layout(
main.title = "Population Growth of States (2001-2011)",
main.title.position = c("center"),
main.title.size = 1,
legend.position = c("right", "bottom")
) +
tm_credits("Data:\n2001-2011 Census", position = c("left", "bottom"))
#to see the map called growth just created above
growth## Variable "decadal_growth" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
#to create a density map
density <- tm_shape(states_sf) + # write the name of the map file
tm_fill(col = "density_km2", title = "No. People / Sq Km",
palette = "YlGnBu") +
tm_borders(lwd = 0.5) +
tm_layout(
main.title = "Population Density of States (2011)",
main.title.position = c("center"),
main.title.size = 1,
legend.position = c("right", "bottom")
) +
tm_credits("Data:\n2011 Census", position = c("left", "bottom"))#to see both the maps created together
tmap_arrange(growth, density)## Variable "decadal_growth" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable "decadal_growth" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
# Inset map creation -------------------------------------------------------
#three maps are to be created..map of india, north east and sex ratio
#to change sequence of palette colour, use a minus "-" sign
#mutate command is used to make new variable..similar to transform variable
#NorthEast Sex and inset map
ne_sex <- simp_sf %>%
filter(region == "Northeastern") %>%
tm_shape() +
tm_fill(col = "sex_ratio", title = "Sex Ratio", palette = "-Reds") +
tm_borders(lwd = 0.5) +
tm_text('state_ut', size = 0.75) +
tm_layout(
main.title = "Sex Ratio in India's Northeast",
main.title.position = c("center"),
main.title.size = 1
) +
tm_credits("Data Source: Wikipedia", position = c("left", "top"))
# To create regional map
regional_sf <- simp_sf %>%
group_by(region) %>%
summarise(pop = sum(pop_2011))
head(regional_sf) ## Simple feature collection with 6 features and 2 fields
## geometry type: GEOMETRY
## dimension: XY
## bbox: xmin: 71.88402 ymin: 6.780492 xmax: 97.4091 ymax: 37.04546
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 3
## region pop geometry
## <chr> <dbl> <GEOMETRY [°]>
## 1 Arabian Sea 6.44e4 MULTIPOLYGON (((71.91522 12.38769, 71.92342 12.34714~
## 2 Bay of Ben~ 3.80e5 MULTIPOLYGON (((92.87543 12.88287, 92.9266 12.88573,~
## 3 Central 9.81e7 POLYGON ((82.81637 23.95908, 82.95171 23.87095, 83.0~
## 4 Eastern 2.70e8 MULTIPOLYGON (((88.88482 21.52353, 88.8334 21.5402, ~
## 5 Northeaste~ 4.56e7 MULTIPOLYGON (((88.75514 27.14671, 88.70399 27.15326~
## 6 Northern 3.00e8 POLYGON ((78.9706 31.23303, 78.9749 31.30516, 79.013~# create inset map
inset <- regional_sf %>%
filter(!region == "Arabian Sea",
!region == "Bay of Bengal") %>%
mutate(northeast = ifelse(region == "Northeastern", TRUE, FALSE)) %>%
tm_shape() +
tm_fill(col = "northeast", palette = c("grey", "red")) +
tm_style("cobalt") +
tm_legend(show = FALSE)# print map using grid
library(grid)
ne_sex
print(inset, vp = viewport(0.24, 0.18, width = 0.2, height = 0.4))
# Creation of faceted maps
# create custom labels for log scale
gdp_seq <- 10 ^ (seq(2.8, 4.0, by = 0.2))
gdp_vec <- scales::dollar(round(gdp_seq))
my_labels = vector(mode = "character", length = 6)
for (i in seq_along(1:6)) {
my_labels[i] = str_c(gdp_vec[i], " to ", gdp_vec[i + 1])
}
simp_sf %>%
mutate(
log_pc_usd = log10(per_capita_gdp_usd),
region_fac = factor(region, levels = c("Northern", "Western", "Southern",
"Central", "Eastern", "Northeastern",
"Arabian Sea", "Bay of Bengal"))
) %>%
filter(!state_ut %in% c("Andaman & Nicobar Islands",
"Lakshadweep")) %>%
tm_shape() +
tm_borders(lwd = 0.5, col = "white") +
tm_fill(col = 'log_pc_usd', title = '', palette = "viridis",
labels = my_labels) +
tm_facets(by = "region_fac", nrow = 2, free.coords = TRUE) +
tm_layout(
main.title = "Per Capita GDP by Region",
main.title.size = 1,
main.title.position = "center",
legend.outside.position = "right"
) +
tm_compass(type = "8star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), size = 1)
# Creation of proportional maps
pop_bubbles <- simp_sf %>%
tm_shape() +
tm_polygons() +
tm_bubbles(col = "gold", size = "pop_2011",
scale = 3, title.size = "") +
tm_text("abb", size = "pop_2011", root = 5,
legend.size.show = FALSE) +
tm_layout(
main.title = "Population (2011)",
main.title.position = c("center"),
main.title.size = 1,
legend.position = c("right", "bottom")
) +
tm_compass(type = "8star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), size = 1)
pop_bubbles## Scale bar set for latitude km and will be different at the top and bottom of the map.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.6. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## The legend is too narrow to place all symbol sizes.
gdp_bubbles <- simp_sf %>%
tm_shape() +
tm_polygons() +
tm_bubbles(col = "gold", size = "nominal_gdp_usd",
scale = 3, title.size = "") +
tm_text("abb", size = "nominal_gdp_usd", root = 5,
legend.size.show = FALSE) +
tm_layout(
main.title = "Nominal GDP (USD)",
main.title.position = c("center"),
main.title.size = 1,
legend.position = c("right", "bottom")
) +
tm_compass(type = "8star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), size = 1)
gdp_bubbles## Scale bar set for latitude km and will be different at the top and bottom of the map.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.52. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## The legend is too narrow to place all symbol sizes.
tmap_arrange(pop_bubbles, gdp_bubbles)## Scale bar set for latitude km and will be different at the top and bottom of the map.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.49. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## The legend is too narrow to place all symbol sizes.## Scale bar set for latitude km and will be different at the top and bottom of the map.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.42. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## The legend is too narrow to place all symbol sizes.
#### Using GGPLOT for maps
library(ggplot2)
library(ggrepel)
st_crs(simp_sf)## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"proj_sf <- simp_sf %>%
st_transform(crs = 24343) %>%
mutate(
CENTROID = purrr::map(geometry, st_centroid),
COORDS = purrr::map(CENTROID, st_coordinates),
COORDS_X = purrr::map_dbl(COORDS, 1),
COORDS_Y = purrr::map_dbl(COORDS, 2)
)
kerala <- proj_sf %>%
filter(state_ut == "Kerala")
head(kerala)## Simple feature collection with 1 feature and 21 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 485275.7 ymin: 916462.7 xmax: 764601.8 ymax: 1412551
## epsg (SRID): 24343
## proj4string: +proj=utm +zone=43 +a=6377299.151 +b=6356098.145120132 +towgs84=295,736,257,0,0,0,0 +units=m +no_defs
## state_ut abb type data_year comparable_economy nominal_gdp_usd
## 1 Kerala KL State 2018–19 Finland 1.2e+11
## nominal_gdp_inr pop_2011 decadal_growth rural_pop urban_pop area_km2
## 1 7.73e+12 33387677 0.049 17445506 15932171 38863
## density_km2 sex_ratio region per_capita_gdp_inr per_capita_gdp_usd
## 1 859 1084 Southern 231522.5 3594.14
## CENTROID COORDS COORDS_X COORDS_Y
## 1 654330.2, 1154519.1 654330.2, 1154519.1 654330.2 1154519
## geometry
## 1 MULTIPOLYGON (((656577.1 12...plot(kerala)## Warning: plotting the first 10 out of 21 attributes; use max.plot = 21 to
## plot all
proj_sf %>%
filter(!state_ut %in% c("Daman & Diu", "Dadra & Nagar Haveli")) %>%
ggplot() +
geom_sf(aes(fill = sex_ratio), lwd = 0) +
geom_sf(fill = NA, color = "grey", lwd = 0.5) +
scale_fill_viridis_c("Sex Ratio", labels = scales::comma, option = "A") +
labs(
title = "Sex Ratio across Indian States",
caption = "Source: Wikipedia"
) +
geom_text_repel(
data = kerala,
mapping = aes(x = COORDS_X, y = COORDS_Y, label = state_ut),
nudge_x = -0.5,
nudge_y = -1
) +
scale_y_continuous(NULL) +
scale_x_continuous(NULL) +
theme(plot.title = element_text(hjust = 0.5)) +
# remove graticules
coord_sf(datum = NA) +
theme_void()
### Dot density maps
# save geometry
proj_geometry <- proj_sf %>% select(state_ut)
#plot(proj_geometry)
# gather data and rejoin geometry
pop_gathered <- proj_sf %>%
st_set_geometry(NULL) %>%
select(state_ut, rural_pop, urban_pop) %>%
gather(key = "pop", value = "count", -state_ut) %>%
arrange(state_ut) %>%
left_join(proj_geometry) %>%
st_as_sf()## Joining, by = "state_ut"# create a list of urban and rural populations
pop_split <- pop_gathered %>% split(.$pop)
# draw 1 dot per 1 lakh people
generate_samples <- function(data) {
st_sample(data, size = round(data$count / 1e5))
}
# generate samples for each and combine
points <- map(pop_split, generate_samples)
points <- imap(points, ~st_sf(tibble(
pop = rep(.y, length(.x))),
geometry = .x))
points <- do.call(rbind, points)
# group points into multipoints
points <- points %>%
group_by(pop) %>%
summarise()# plot with ggplot
points %>%
ggplot() +
geom_sf(data = simp_sf) +
geom_sf(aes(color = pop, fill = pop),
size = 0.1, alpha = 0.4) +
scale_fill_discrete("Population", labels = c("Rural", "Urban")) +
labs(
title = "Density of India's Urban and Rural Population (2011)",
caption = "1 dot = 1 lakh people"
) +
theme(plot.title = element_text(hjust = 0.5)) +
coord_sf(datum = NA) +
theme_void() +
guides(color = FALSE)
# End ---------------------------------------------------------------------#3 types of spatial data #geostatistical, regional and point data #fallacies..ecological and Modifiable Unit Area Problem (MUAP)
#to save file in ESRI #check and confirm that class is sf #proj_sf #class(proj_sf) #to save as shp file #st_write(proj_sf, “proj_sf_shp.shp”)