google plus level 2 polygons extracted from file previously compiled for MENA/global study

mk_gpl2_points_ETH.R

# google plus level 2 polygons extracted from file previously compiled for MENA/global study
# source https://grid.plus.codes/

setwd("D:/sm2/MENA/ge")

reg <- "AFR"
iso = "ETH"

x <- as.list(c("sp","sf","rgdal","rgeos","raster","readxl","dplyr","data.table"))
lapply(x, require, character.only = TRUE)

#get WB_Name and Region
gaul_wb <- read_excel("gaul_adm0.xls",sheet="gaul_adm0")
gaul_wb_sel <- gaul_wb %>%
  dplyr::select(-Names,-First_ADM0,-Regions,-hs_region) %>%
  rename(WB_Region = hs_regco,WB_ISO = ISO_Codes) %>%
  filter(WB_ISO == iso)

#get admin boundaries
adm <- sf::st_read(dsn = "D:/sm2/MENA/ge",
                   layer = "g2015_2014_2",stringsAsFactors=F)
adm <- adm[,c("ADM2_CODE","ADM1_CODE","ADM0_CODE")]
adm <- merge(adm,gaul_wb_sel,all.x=F) #drops countries outside region

#get gp level2 polygons for reg
plys_eh <- st_read("gp2_EH.shp")
#save(plys_eh,file="reg_rdat/gp2_EH.RData")
load("reg_rdat/gp2_EH.RData")
plys <- plys_eh[adm,]
st_write(plys,paste0("gpl2_",iso,".shp")

#make centerpoints for sampling at level2 (0.05 deg) resolution
gpl2_ctr <- st_centroid(plys)
ctradm <- st_join(gpl2_ctr,adm, join=st_intersects, left=F) 
ctradm$lon <- st_coordinates(ctradm)[1]
ctradm$lat <- st_coordinates(ctradm)[2]
ctradm_df <- as.data.frame(st_set_geometry(ctradm, NULL))
 
#then loop thru cells make points at 0.01 deg spacing, 25 per cell
#and assign ctr, gp2 name

oldw <- getOption("warn")
options(warn = -1)

samples <- list()
  for (i in 1:nrow(plys)) {
    psub <- plys[i,]
    lon <- seq(extent(psub)[1]-0.005,extent(psub)[2],0.01)
    lon <- lon[-1]
    lat <- seq(extent(psub)[3]-0.005,extent(psub)[4],0.01)
    lat <- lat[-1]
    pts <- expand.grid(lon,lat)
    coordinates(pts) <- ~Var1+Var2
    crs(pts) <- crs(adm)
    samples[[i]] <- sp::SpatialPointsDataFrame(pts, 
      data.frame(Name = rep(psub$Name,length(pts))))
#                 WB_ISO = rep(psub$WB_ISO,length(pts))))
    if (i %% 1000==0){
      print(paste0("completed ",i," at ",Sys.time()))
    }
    flush.console()
  }

options(warn = oldw)

dfouts <- list()
for(i in 1:length(samples)) { dfouts[[i]] <- as.data.frame(samples[[i]]) }
ptsall <- rbindlist(dfouts)
names(ptsall)[2] <- "lon"
names(ptsall)[3] <- "lat"
ptsall$pid <- seq_len(nrow(ptsall))
ptsall <- ptsall[,c(4,1,2,3)]

#join adm
ptsall_sf <- st_as_sf(ptsall, coords = c("lon", "lat"), crs = crs(adm))
ptsadm <- st_join(ptsall_sf,adm, join=st_intersects, left=F)
ptsadm$lon <- st_coordinates(ptsadm)[1]
ptsadm$lat <- st_coordinates(ptsadm)[2]
ptsadm_df <- as.data.frame(st_set_geometry(ptsadm, NULL))

#write outputs
fwrite(ctradm_df,paste0("shp_out/gp_level2_ctr_",iso,".csv"))
st_write(ctradm, paste0("shp_out/gpl2_ctr_",iso,".shp"))
fwrite(ptsadm_df,paste0("shp_out/gp_level2_25pts_",iso,".csv"))
st_write(ptsadm, paste0("shp_out/gpl2_25pts_",iso,".shp"))