Generate SpatialPointDataFrame based on dataset in Xitou

xitou_gis.R

#!/usr/bin/env Rscript
#
# ----- For NTU Forestry 3049 Course Only -----
#
# Generate SpatialPointDataFrame based on dataset in Xitou
#
# This Gist: https://git.io/v6ozP
#
# Mail:     b01605002@ntu.edu.tw
# Author:   Weihang Lo
# Date:     2016.08


library(data.table)
library(readxl)
library(rgdal)
library(rgeos)


#---------------------------------------
# Helper Functions
#---------------------------------------

# general function: affine transformation
affineTransform <- function(v11, v12, v21, v22, v31, v32, pts) {
    transMat <- matrix(c(v11, v12, 0, v21, v22, 0, v31, v32, 1), ncol = 3)
    (transMat %*% rbind(pts, 1))[1:2, , drop = FALSE]
}

# general function: coordinates transformation (using affine transform)
coordRotate <- function(coords, origin, angle) {
    ang <- angle / 180
    coords_t <- t(coords)
    if (dim(coords_t)[1] == 1)
        coords_t <- t(coords_t)
    translated <- affineTransform(1, 0, 0, 1, -origin[1], -origin[2], coords_t)
    t(affineTransform(cospi(ang), -sinpi(ang), sinpi(ang), cospi(ang),
                      origin[1], origin[2], translated))
}


# helper function: cleaning dirty excel data
cleanDataFrameToDataTable <- function(x, na.rm = FALSE) {
    if (na.rm)
        x <- x[complete.cases(x), ]
    setDT(x)
    setnames(x, sapply(names(x), function(x) gsub("\\s|\\.", "_", x)))
    return(x)
}

# helper function: enlarge bounds of spatial object for rotating grid
generateCirCularBounds <- function(spgeom) {
    midPt <- rowMeans(spgeom@bbox)
    midPt <- SpatialPoints(matrix(midPt, ncol = 2),
                           proj4string = spgeom@proj4string)
    gBuffer(midPt, width = dist(t(spgeom@bbox))/ 2, quadsegs = 256)
}

# helper function: sample valid points in spatial polygon
samplePts <- function(spgeom, ntree, type, iter = 20) {
    if (type == 'natural') {
        pts <- spsample(spgeom, ntree, type = 'random', iter = iter)
        pts@proj4string <- spgeom@proj4string
        return(pts)
    }

    bounds <- generateCirCularBounds(spgeom)
    sampleSize <- max(gArea(bounds) / gArea(spgeom) * ntree, 1)

    i <- 0
    result <- NULL
    resultDelta <- NULL
    while (i < iter) {
        pts <- spsample(bounds@polygons[[1]], n = sampleSize,
                        type = type, iter = iter)
        if(is.null(pts))
            next

        pts@proj4string <- spgeom@proj4string

        midPt <- rowMeans(spgeom@bbox)
        pts@coords <- coordRotate(pts@coords, midPt, runif(1) * 359)
        pts <- pts[rowSums(is.na(pts %over% spgeom)) == 0, ]
        delta <- abs(length(pts) - ntree)

        if (is.null(result) || delta < resultDelta) {
            result <- pts
            resultDelta <- delta
        }

        if (resultDelta < min(10, ntree * 0.95))
            break

        i = i + 1
    }

    return(result)
}


# helper funciton: sampling dbh with normal distribution
sampleDBH <- function(f, s, size, data) {
    dataDBH <- data[forest == f & stand == s,
                        .(dbh_mean, dbh_sd, dbh_min, dbh_max)]
    newDBHs <- dataDBH[, .(rnorm(size, dbh_mean, dbh_sd))]
    setnames(newDBHs, "dbh")
    newDBHs[dbh < dataDBH$dbh_min, dbh := dataDBH$dbh_min]
    newDBHs[dbh > dataDBH$dbh_max, dbh := dataDBH$dbh_max]
    return(newDBHs)
}

# helper function: final aggregate function for generate trees
# param spgeom:     spatial object
# param data:       forest/stand data
# param type:       natural or plantation
# param compData:   species composition data
generateTrees <- function(spgeom, data, type, compData) {
    coordList <- vector(mode = "list", nrow(data))
    for (i in seq_len(nrow(data))) {
        cat("index ", i, "\n")

        isNatural <- type == "natural"

        f <- data[i, forest]
        s <- data[i, stand]

        randPts <- samplePts(spgeom[i, ], data[i, ntree],
                             ifelse(isNatural, "random", "regular"))
        spdf <- data.table(x = randPts@coords[, 1], y = randPts@coords[, 2])
        spdf[, c("forest", "stand") := list(f, s)]
        spdf[, dbh := sampleDBH(f, s, length(randPts), data)]
        if (isNatural) {
            comp <- compData[forest == f & stand == s, .(species, proportion)]
            spp <- sample(comp$species, length(randPts), TRUE, comp$proportion)
            spdf[, species := .(spp)]
        } else {
            spp <- data[i, species]
            spdf[, species := .(data[i, species])]
        }

        coordList[[i]] <- spdf
    }
    coords <- rbindlist(coordList)
    trees <- SpatialPointsDataFrame(coords[, .(x, y)],
                data = coords[, -c("x", "y"), with = FALSE])
    trees@data <- cleanDataFrameToDataTable(trees@data, na.rm = TRUE)
    return(trees)
}


#---------------------------------------
# Handle Real Dataset
#---------------------------------------

# Calculate area per polygon in hectare (ha)
xitou <- readOGR(dsn = "Data/GIS/", layer = "xitou")
setDT(xitou@data)
xitou@data$area <- sapply(xitou@polygons, function(x) x@area / 1e4)

# Set Coordinate Reference System
TWD97_TM2_ZONE121 <- CRS("+proj=tmerc 
                         +lat_0=0 +lon_0=121 
                         +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 
                         +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
xitou@proj4string <- TWD97_TM2_ZONE121

# Natural Forest -----------------------
naturalExcel <- "Data/Xitou/xitou_natural.xlsx"

naturalData <- read_excel(naturalExcel, "base data")
naturalData <- cleanDataFrameToDataTable(naturalData, na.rm = TRUE)
naturalComp <- read_excel(naturalExcel, "species composition")
naturalComp <- cleanDataFrameToDataTable(naturalComp, na.rm = TRUE)

# join spatial data
naturalData <- merge(naturalData, xitou, by = c("forest", "stand"))
naturalData[, ntree := tree_ha * area]

# generate points in polygons
naturalTrees <- generateTrees(xitou, naturalData, 'natural', naturalComp)
naturalTrees@proj4string <- TWD97_TM2_ZONE121
writeOGR(naturalTrees, ".", "naturalTrees", driver = "ESRI Shapefile")


# Plantation Forest --------------------

plantationExcel <- "Data/Xitou/xitou_plantation.xlsx"

plantationData <- read_excel(plantationExcel, "base data")
plantationData <- cleanDataFrameToDataTable(plantationData, na.rm = TRUE)
plantationData[is.na(plantationData)]  <- 1.1

# join spatial data
plantationData <- merge(plantationData, xitou, by = c("forest", "stand"))
plantationData[, ntree := tree_ha * area]

# generate points in polygons
plantationTrees <- generateTrees(xitou, plantationData, 'plantation')
plantationTrees@proj4string <- TWD97_TM2_ZONE121
writeOGR(plantationTrees, ".", "plantationTrees", driver = "ESRI Shapefile")