Spatial data analysis and visualization play a crucial role in understanding the geographic landscape of a region. In this tutorial, we’ll see how to create detailed maps of Bangladesh at the division and district levels using the sf (Simple Features) package in R.

The sf package offers a comprehensive set of tools for working with spatial data, making it an excellent choice. The first is to load the necessary libraries and download the shapefile (I used this website: https://www.diva-gis.org/gdata)

# read the shapefile at division level
sf1 <- st_read(dsn="~/Downloads/BGD_adm/BGD_adm1.shp", layer="BGD_adm1")
## Reading layer `BGD_adm1' from data source 
##   `/Users/macintosh/Downloads/BGD_adm/BGD_adm1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 88.01057 ymin: 20.74111 xmax: 92.67366 ymax: 26.6344
## Geodetic CRS:  WGS 84
# read the shapefile at district level
sf2 <- st_read(dsn="~/Downloads/BGD_adm/BGD_adm2.shp", layer="BGD_adm2")
## Reading layer `BGD_adm2' from data source 
##   `/Users/macintosh/Downloads/BGD_adm/BGD_adm2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 65 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 88.01057 ymin: 20.74111 xmax: 92.67366 ymax: 26.6344
## Geodetic CRS:  WGS 84

Right now, there is no real-life variable to plot in the data. So, we will make a fictitious population column to make the plot.

First, we will make the plot at the division level:

sf1$population <- sample(50000:500000, nrow(sf1), replace = T)

ggplot() + 
    geom_sf(data = sf1, aes(fill = population)) + 
    scale_fill_gradient(low = "#f8ead6", high = "#98df8a")


# And then at the district level:

sf2$population <- sample(50000:500000, nrow(sf2), replace = T)

ggplot() + 
    geom_sf(data = sf2, aes(fill = population)) + 
    scale_fill_gradient(low = "#f8ead6", high = "#98df8a")