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.


#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.


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.



#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.

  • Thanks! I want to keep the data in a grid of the same level of resolution (not aggregated) exploitable for statistical analysis (with spml or spdep I guess although I haven't got that far yet). I guess I need the data in vector format for this, not raster format. Sorry if my question isn't clear as I'm a beginner at this.
    – syre
    May 28 '20 at 3:32
  • sca_gadm is already in the same CRS (WGS84) as the raster so I guess there is no need for reprojection in your approach? However, I think I need planar coordinates for later stages.
    – syre
    May 28 '20 at 3:43
  • 1
    Please note that the velox package is depreciated. May 28 '20 at 5:38
  • 2
    An alternative to velox is the exactextractr package. It is actually faster than velox and coercing to a different object class is not necessary as the exact_extract function operates on raster class objects. Also, since coercion is not necessary for exact_extract, cropping the raster data to the extent of the polygons is not necessary and would actually slow things down. Jun 1 '20 at 16:04
  • 1
    It will be notability faster to project the polygons to the same CSR as the raster, defining a new polygon object and then join the tabular results back to the original polygons. Although, it is good practice to define a common projection and get all of your data correctly aligned to a standard analysis extend. This is part of good data management and ensures reproducibility down the line. Jun 2 '20 at 14:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.