R/gridPoly_OverlapFraction.R

In kongdd/CMIP5tools: Functions to Tidy and Analysis CMIP5/6 Data

#' faction of grids of SpatialPixels located in polygon
#' 
#' @param grid SpatialGrid, with the column of `id`.
#' @param poly SpatialPolygon obj 
#' @param write.fig boolean
#' 
#' @return
#' * `grid`    : original `grid`
#' * `gridclip`: filter by fraction of clip area
#' * `I_sel`   : selected gridId
#' * `Id`      : Id of original `grid`
#' 
#' @examples
#' \dontrun{
#' gridInfo <- gridPoly_OverlapFraction(gridclip, poly_CH)
#' d <- gridInfo$grid@data %>% data.table()
#' }
#' 
#' @export
gridPoly_OverlapFraction <- function(grid, poly, write.fig = FALSE){
    Id0 <- grid$id
    ## why st not work?
    # st_poly <- st_as_sfc(poly)
    # st_poly_grid <- st_as_sfc(poly_grid)
    grid$I    <- 1:nrow(grid)
    poly_grid <- as(grid, "SpatialPolygonsDataFrame")
    poly_clip <- raster::intersect(poly_grid, poly)
    
    sp_layout <- list("sp.polygons", poly, fill = "transparent", 
        lwd = 0.5, first = FALSE)
    # sp_poly <- list("sp.polygons", poly, first = FALSE) 
    
    I <- poly_clip$I #%>% unique()
    # union polygons
    AREA_clip <- raster::area(poly_clip) # not aggregated
    I_overlap <- unique(I) # combine duplicate
    I_sel0    <- I_overlap # select all overlap grids

    m2_to_km2 <- 1/1e6
    s_clip    <- split(AREA_clip, I) %>% map_dbl(sum, na.rm = TRUE)    
    s_grid <- raster::area(poly_grid[I_overlap, ])
    fraction  <- pmin(s_clip/s_grid, 1)
    poly_clip$fraction <- AREA_clip / raster::area(poly_grid[I, ])
    
    # I_sel0 <- I_overlap[fraction > value]
    ## 1. return fraction to grid
    area_0 <- fraction_0 <- rep(NA, length(grid))
    w <- rep(0, length(grid))
    
    area_0[I_overlap]     <- s_clip   # note that: it is clip_region area, not grid area
    fraction_0[I_overlap] <- fraction
    # w            <- area_0 / sum(area_0, na.rm = TRUE) # weights of every pixel
    w[I_sel0]      <- (area_0/sum(area_0, na.rm = TRUE))[I_sel0] ## weights by area
    # w[!is.finite(w)] <- 0

    ## select all overlap grids
    grid$fraction <- fraction_0
    grid$area     <- area_0 * m2_to_km2
    grid$w        <- w
    
    ## 2. Compare different fraction
    # spplot(poly_clip, "fraction",
    #    col.regions = RColorBrewer::brewer.pal(9, "Blues"),
    #    # main = title,
    #    sp.layout = sp_layout)
    
    brks <- seq(0, 1, 0.1)
    blues = c("#F7FBFF", "#DEEBF7", "#C6DBEF", "#9ECAE1", "#6BAED6", "#4292C6", 
        "#2171B5", "#08519C", "#08306B")
    colfun <- colorRampPalette(blues) # RColorBrewer::brewer.pal(9, "Blues") 
    colors <- colfun(length(brks)-1)

    range <- grid@bbox %>% t() %>% as.numeric()
    xlim  <- range[1:2]
    ylim  <- range[3:4]

    grps = c(0, 0.1, 0.2, 0.3)[1]
    ps <- foreach(value = grps) %do% {
        title <- sprintf("fraction >= %d%%", value*100)
        I_sel <- I_overlap[fraction > value]
        # pi <- poly[, ]
        spplot(grid[I_sel, ], "fraction", 
               col.regions = colors, at = brks,
               main = title, xlim = xlim, ylim = ylim,
               lwd = 0.5, col = "black",
               sp.layout = sp_layout)
    }
    
    if (write.fig) {
        g <- gridExtra::grid.arrange(grobs = ps, nrow = 2) # arrangeGrob
        # write_fig(g, "Figure/Select_grid_v2.pdf", 10, 7)
    }
    
    ## 3. select
    grid@data %<>% data.table()
    gridclip <- grid[I_sel0, ]
    # listk(grid, gridclip, I_sel = I_sel0)
    gridclip
}