## Quadtree-based decomposition of rasters (gridded population dataset)
## Required packages
install.packages("raster")
install.packages("maptools")
install.packages("rgdal")
install.packages("sp")
install.packages("doSNOW")
## Set directory
setwd("C:/")
## Define the raster data (Gridded population data)
RasataData=raster("Raster data")
## Define the shape file to identfy the adminstrative regions or settelment types.
reg=readOGR("shape file")
RasterData.Region1=subset(reg, ID=="First Region name")
RasterData.Region1=subset(reg, grepl("Second Region name", ID))
Awdal.rast=crop(som, som.Awdal)
Banaadir.rast=crop(som, som.Banaadir)
## Regions to process
provs <- as.character(reg$ID)
## help function: split images into quads
splitim <- function(im){
r <- nrow(im)
c <- ncol(im)
dr <- round(r/2)
dc <- round(c/2)
quads <- list()
quads[[1]] <- im[1:dr, 1:dc, drop=FALSE]
quads[[2]] <- im[(dr+1):r, 1:dc, drop=FALSE]
quads[[3]] <- im[(dr+1):r, (dc+1):c, drop=FALSE]
quads[[4]] <- im[1:dr, (dc+1):c, drop=FALSE]
return(quads)
}
## Pre-define parameters for population size and area
minSize <- 3500 # Maximum population size
gridSize <- 32 # maximum grid dimension
### main loop
for(p in provs[1:18]){
print(p)
grid_c<- crop(som, subset(reg, ID==p))
grid<-mask(grid_c, subset(reg, ID==p))
#grid=mask(grid, UR)
dim(grid)
iml <- list() # store the list of images to process
iml[[1]] <- grid # first in the image list is the full grid
res <- FALSE # results, start false to begin with a split
res2 <- FALSE
imS <- list() # output list of split grids
# main processing loop
while (prod(res)==0 | prod(res2)==0) {
iml <- unlist(lapply(iml, splitim))
res <- unlist(lapply(iml, function(i){ cellStats(i, stat='sum') < minSize } ))
res2 <- unlist(lapply(iml, function(i){ max(dim(i)) < gridSize } ))
imS <- c(imS, iml[which(res==1 & res2==1)])
iml <- iml[which(res==0 | res2==0)]
}
## drop subscenes with no data
keep <- lapply(1:length(imS), function(i){ cellStats(imS[[i]],'mean')>0 & !is.nan(cellStats(imS[[i]], 'mean')) })
imS <- imS[unlist(keep)==T]
## output split grids
if(!dir.exists(paste0( "subscenes"))){dir.create(paste0( "subscenes"))}
if(!dir.exists(paste0( "subscenes/", p,"/"))){dir.create(paste0( "subscenes/", p,"/"))}
## output extent polygons
plylist <- lapply(1:length(imS),
function(i){
if(!is.null(imS[[i]])){
pop=cellStats(imS[[i]], sum)
e <- extent(imS[[i]])
sp <- as(e, 'SpatialPolygons')
crs(sp) <- crs(imS[[i]])
polydf <- SpatialPolygonsDataFrame(sp, data=data.frame(id=i, pop=pop))
return(polydf)}
})
## merge all to a single feature
joinpoly <- do.call(bind, plylist)
if(!dir.exists(paste0( "subshapes"))){dir.create(paste0( "subshapes"))}
if(!dir.exists(paste0( "subshapes/", p,"/"))){dir.create(paste0( "subshapes/", p,"/"))}
shapefile(joinpoly, paste0( "subshapes/", p, "/", p, "_ply_all3.shp"), overwrite=T)
pdf(paste0( "subshapes/", p, "/", p, "_plot.pdf"))
plot(grid, main=p)
plot(joinpoly, add=T)
plot(reg, add=T)
dev.off()
}
Qader, S.H., Lefebvre, V., Tatem, A.J. et al. Using gridded population and quadtree sampling units to support survey sample design in low-income settings. Int J Health Geogr 19, 10 (2020). https://doi.org/10.1186/s12942-020-00205-5
.
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent