library(tidyverse)
library(data.table)
library(sf)
library(tmap)
library(readr)
library(geojsonsf)
library(osmdata)
library(basemapR)
library(RColorBrewer)
library(classInt)
library(R.utils)
library(dplyr)
library(ggplot2)
library(viridis)
library(raster)
library(terra)
library(exactextractr)
library(tidyterra)
library(classInt)Assignment_1
Installing Packages
Your Annotated Code
raster_2010 <- rast("/Users/Maciej/Desktop/GDS/gpw-v4-population-count-rev11_2010_2pt5_min_tif/gpw_v4_population_count_rev11_2010_2pt5_min.tif") # import raster data 2010
raster_2020 <- rast("/Users/Maciej/Desktop/GDS/gpw-v4-population-count-rev11_2020_2pt5_min_tif/gpw_v4_population_count_rev11_2020_2pt5_min.tif") # import raster data 2020
china_shp <- read_sf("/Users/Maciej/Desktop/GDS/gadm36_CHN_shp/gadm36_CHN_2.shp") # import china shapefile
plot(china_shp$geometry) # plot china boundary
crs(china_shp) # check CRS assigned to china_shp[1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"crs(raster_2010) # check CRS assigned to raster_2010[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"crs(raster_2020) # check CRS assigned to raster_2020[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"current_crs <- st_crs(china_shp)
new_crs <- st_crs("+init=epsg:4490") # define the new CRS EPSG:4490
china_shp <- st_set_crs(china_shp, new_crs) # change the CRS of 'china_shp' to the new CRS
new_crs <- st_crs(china_shp) # check the new CRS
raster_2010 <- terra::project(raster_2010, crs(china_shp)) # reporjectig raster_2010 data to the crs of the china_shp
|---------|---------|---------|---------|
=========================================
crs(raster_2010)[1] "GEOGCRS[\"China Geodetic Coordinate System 2000\",\n DATUM[\"China 2000\",\n ELLIPSOID[\"CGCS2000\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",1043]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]],\n CS[ellipsoidal,2],\n AXIS[\"longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"China - onshore and offshore.\"],\n BBOX[16.7,73.62,53.56,134.77]]]"raster_2020 <- terra::project(raster_2020, crs(china_shp)) # reporjectig raster_2020 data to the crs of the china_shp
|---------|---------|---------|---------|
=========================================
crs(raster_2020)[1] "GEOGCRS[\"China Geodetic Coordinate System 2000\",\n DATUM[\"China 2000\",\n ELLIPSOID[\"CGCS2000\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",1043]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]],\n CS[ellipsoidal,2],\n AXIS[\"longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"China - onshore and offshore.\"],\n BBOX[16.7,73.62,53.56,134.77]]]"pop_2010 <- crop(raster_2010, extent(china_shp)) #changes extent of the raster
pop_2010_mask <- mask(pop_2010, china_shp) # removes other countries
pop_2020 <- crop(raster_2020, extent(china_shp)) #changes extent of the raster
pop_2020_mask <- mask(pop_2020, china_shp) # removes other countries
china_shp_area <- st_area(china_shp)*0.000001 # converts square meters to square kilometers for china boundary
population_2010 <- raster::extract(pop_2010_mask, china_shp, fun = sum, na.rm = T) #extracts data from each geometry in china_shp
population_2020 <- raster::extract(pop_2020_mask, china_shp, fun = sum, na.rm = T) #extracts data from each geometry in china_shp into new variable
head(population_2020) ID gpw_v4_population_count_rev11_2020_2pt5_min
1 1 5509757
2 2 3048208
3 3 4665156
4 4 4084265
5 5 1420530
6 6 4035317pop_density_2020 <- population_2020$gpw_v4_population_count_rev11_2020_2pt5_min / china_shp_area # perform a calculation to compute population density for the year 2020
head(pop_density_2020)Units: [1/m^2]
[1] 357.2255 516.0824 542.7450 432.4352 167.6626 297.2462china_shp$Density_per_sq_km <- pop_density_2020 # assign the values of the variable pop_density_2020 to a new column called Density_per_sq_km
summary(china_shp$Density_per_sq_km) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.416 122.207 293.594 439.135 564.267 6476.750 tm_basemap() +
tm_shape(china_shp) +
tm_polygons("Density_per_sq_km", palette = "YlGn", id = "VARNAME_2", n = 5, style = "quantile") +
tm_scale_bar(position = c("right", "bottom"), text.size = 0.5, color.dark = "black", color.light = "white") +
tm_compass(type = "arrow", position = c("left", "top"), size = 1.5) +
tm_layout(main.title = "China population density 2020", main.title.position = c("center", "top"), legend.frame = TRUE, legend.height = 0.25)
## Population change
china_pop_2010 <- raster(pop_2010_mask) # create a new raster object named china_pop_2010 by creating raster based on pop_2010_mask
china_pop_2020 <- raster(pop_2020_mask)
population_change <- overlay(china_pop_2010, china_pop_2020, fun = function(x,y)y-x) # calculate population change between china_pop_2010 and china_pop_2020
summary(population_change) # get statistick for population_change layer
Min. -16164.593750
1st Qu. 0.000000
Median 2.639648
3rd Qu. 51.269241
Max. 217922.500000
NA's 694030.000000reclass_rules <- matrix(c(-Inf, 0, 0, 2.645633, 2.64533, 51.306091, 51.306091, Inf, 1, 2, 3, 4), ncol = 2, byrow = TRUE) # reclassify data values using min, mean, 1st Qu, median, 3rd Qu, max
reclassified_raster <- reclassify(population_change, reclass_rules) # perform a reclassification of a raster using predefined rules
tm_shape(reclassified_raster, raster.downsample = FALSE) +
tm_raster(title = "Population change",
palette = c("red", "orange", "lightgreen", "darkgreen"),
breaks = c(0, 1, 2, 3, 4),
labels = c("Decline", "Neutral", "Growth", "High Growth")) +
tm_layout(main.title = "China population change 2010-2020",
main.title.position = c("center", "top"),
legend.position = c("left", "bottom"), legend.frame = TRUE, legend.text.size = 0.7) + # set legend position
tm_scale_bar(position = c("right", "bottom"), text.size = 0.5, color.dark = "black", color.light = "white") +
tm_compass(type = "arrow", position = c("left", "top"), size = 1.5)
# Create the map object by adding layers to the basemap
basemap <- tm_basemap()
map_1 <- basemap +
tm_shape(china_shp) +
tm_polygons("Density_per_sq_km", palette = "YlGn", id = "VARNAME_2", n = 5, style = "quantile") +
tm_scale_bar(position = c("right", "bottom"), text.size = 0.5, color.dark = "black", color.light = "white") +
tm_compass(type = "arrow", position = c("left", "top"), size = 1.5) +
tm_layout(main.title = "China population density 2020", main.title.position = c("center", "top"), legend.frame = TRUE, legend.height = 0.25)
map_2 <- basemap +
tm_shape(reclassified_raster, raster.downsample = FALSE) +
tm_raster(title = "Population change",
palette = c("red", "orange", "lightgreen", "darkgreen"),
breaks = c(0, 1, 2, 3, 4),
labels = c("Decline", "Neutral", "Growth", "High Growth")) +
tm_layout(main.title = "China population change 2010-2020",
main.title.position = c("center", "top"),
legend.position = c("left", "bottom"),
legend.frame = TRUE,
legend.text.size = 0.7) +
tm_scale_bar(position = c("right", "bottom"), text.size = 0.5, color.dark = "black", color.light = "white") +
tm_compass(type = "arrow", position = c("left", "top"), size = 1.5)Maps
map_1
map_2
Up to 500 words which should address the following questions.
- What value is associated with each pixel?..
- Which CRS are you using? Justify your answer.
- What is Map 1 saying about population density in China?
- Which categorical divisions (Decline, Neutral, Growth and High Growth) did you chose for Map 2 and why?
On the China 2020 population density map, every pixel represents a categorical variable allocated to it. These numbers, in this instance, indicate the population density per square kilometre for every multipolygon (district) that is part of the China shapefile. Population changes over time are represented by the values assigned to each pixel, which represent the variations between the raster layer containing data from 2010 and 2020. The CRS used in this project is EPSG 4490. The coordinate reference system is specifically intended for China, which is the primary rationale for utilising these projections instead of the ones that were initially allocated to the data. Compared to WGS 84, which is utilised for worldwide applications, it was achievable to minimise global distortion by using EPSG 4490, which is optimised for accuracy and precision within the geographic boundaries of China. The distribution of population density for every region in China is illustrated in Map 1. The heavily inhabited east and southeast of China have a substantial population density difference from the much less dense west and northwest. The concentration of the most significant agglomerations close to coastal regions undoubtedly influences this distribution of population since it affects the economic development of those areas. Based on the analysis objectives and the interpretation of the population change data, these categories and the corresponding value ranges have been selected. Decline: Population change values falling within the range of (-Inf, 0), Neutral: Population change values in the range of (0, 2.674431), Growth: Population change values in the range of (2.674431, 51.306091), High Growth: Population change values in the range of (51.306091, Inf). These breaks categorise variables into four groups based on quartiles and customised intervals. It makes it easier to identify places with particular characteristics, like decline or rapid growth, by providing a way to express different levels of population change on the map visually.