Convert raster to polygons, crop, reproject

Question

I wish to map a high-resolution grid of population counts of Korea and extract from it the Seoul Capital Area (SCA) at the same level of resolution. The population data comes in two forms: raster (population_kor_2018-10-01_geotiff.zip) or CSV centroid data (population_kor_2018-10-01.csv.zip). The SCA administrative boundaries come from GADM (R (sf) file).

Am I doing this right? in the right order? Is it not correct to reproject the raster before converting to polygons since the variable (population) is areal thus should be recalculated? I have also seen it recommended to reproject polygons rather than the raster.

library(raster)
library(rgdal)
library(sf)

#Load raster of population
kor_dens <- raster("population_kor_2018-10-01.tif")
kor_dens@crs #CRS is WGS84

#Load SCA admin boundaries
sca_gadm <- readRDS("gadm36_KOR_1_sf.rds")
st_crs(sca_gadm)$epsg #CRS is WGS84

#Crop raster down to SCA first to save on computing (rectangular outcome)
kor_dens <- crop(kor_dens, extent(sca_gadm))

#Reproject raster to planar CRS, which adjusts population values
kor_dens <- projectRaster(pop_rast, crs=CRS('+init=EPSG:5179'))

#Reproject SCA boundaries to planar CRS
sca_gadm <- st_transform(sca_gadm, 5179)

#Mask raster, which replaces values outside the area of interest with NA
kor_dens <- mask(kor_dens, as(sca_gadm$geometry, 'Spatial'))

#Convert raster to polygons, where cells with NA are not converted
kor_dens <- rasterToPolygons(kor_dens)

The resulting polygon object is massive at 9.3 Gb, so I'm still not sure this is the way to go for statistical analysis.

EDIT: I've replaced the code entirely, inspired by @BrunoConteLeite's answer.

Answer 1

I would change your approach by using the velox package that does super fast and efficient raster aggregation. I would also avoid re-projecting the raster object as this step is super CPU demanding, but instead, reproject the polygon to the former's projection.

Try:

library(raster)
library(velox)
library(rgdal)
library(sf)

#Load raster
kor_dens <- raster("population_kor_2018-10-01.tif")

sca_gadm <- readRDS("gadm36_KOR_1_sf.rds")
sca_gadm <- sca_gadm[sca_gadm$NAME_1 %in% c("Seoul", "Incheon", "Gyeonggi-do"),]
# make a sp object with which velox works
sca_gadm_sp <- as(sca_gadm, 'Spatial')
# reproject it to the raster's CRS:
sca_gadm_sp <- spTransform(sca_gadm_sp, kor_dens@crs)

# aggregating it:
kor_dens <- crop(kor_dens,extent(sca_gadm_sp)) # cropping the raster over the area of interest
kor_dens <- velox(kor_dens) # making it a velox object
final.df <- kor_dens$extract(sp = sca_gadm_sp, fun = function(x) sum(x, na.rm = TRUE), df=T, small = T) # summing the count of population within the admin boundaries

final.df is the data frame of interest.