vph: Volume Per Hectares Prediction Model

Usage Arguments Examples

Usage

1
vph(data, img, poly, model)

Arguments

data
img
poly
model

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (data, img, poly, model)
{
    require(raster)
    require(rgdal)
    require(rgeos)
    traindata <- base::sample(nrow(data), size = 0.7 * nrow(data),
        replace = FALSE)
    TrainSet <- data[traindata, ]
    ValidSet <- data[-traindata, ]
    if (model == "mod1") {
        mod <- lm(log(VUB) ~ I(log(DBH)) + I(log(H)) + I(log(BA)),
            data = data)
        coef(mod)
        vol <- function(x) {
            exp(coef(mod)[1]) * (23.1^coef(mod)[2]) * (x^coef(mod)[3]) *
                (0.046^coef(mod)[4])
        }
        vol.agg <- aggregate(img, fact = 4/res(img))
        volimg <- raster::calc(vol.agg, vol)
        x <- xres(volimg)
        y <- yres(volimg)
        ha <- (x * y)/10000
        volha.func <- function(x) {
            x/ha
        }
        volha <- raster::calc(volimg, volha.func)
        par(mar = rep(0.5, 4))
        compsRas <- rasterize(poly, volha)
        zoneStat <- zonal(volha, compsRas, "mean")
        poly[["volumeperha"]] <- zoneStat[, "mean"]
        colRamp <- colorRampPalette(c("lightgoldenrod1", "tomato2"))(10)
        polyCols <- colRamp[as.numeric(cut(poly[["volumeperha"]],
            breaks = 10))]
        options(repr.plot.width = 4, repr.plot.height = 4)
        plot(poly, col = polyCols, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n", main = "Volume/Ha")
        text(gCentroid(poly, byid = TRUE), round(poly[["volumeperha"]],
            1), col = "blue", font = 2)
        writeRaster(volha, "./output/volperha_mod1.tif", format = "GTiff",
            overwrite = TRUE)
        writeOGR(poly, "output", "vhamod1", driver = "ESRI Shapefile",
            overwrite = TRUE)
    }
    else if (model == "mod2") {
        mod <- lm(log(VUB) ~ I(log(DBH)) + I(log(H)) + I(log(DBH)^2) +
            I(log(H)^2), data = data)
        coef(mod)
        vol <- function(x) {
            exp(coef(mod)[1]) * (23.1^coef(mod)[2]) * (x^coef(mod)[3]) *
                ((23.1^2)^coef(mod)[4]) * ((x^2)^coef(mod)[5])
        }
        vol.agg <- aggregate(img, fact = 3/res(img))
        volimg <- raster::calc(vol.agg, vol)
        x <- xres(volimg)
        y <- yres(volimg)
        ha <- (x * y)/10000
        volha.func <- function(x) {
            x/ha
        }
        volha <- raster::calc(volimg, volha.func)
        par(mar = rep(0.5, 4))
        compsRas <- rasterize(poly, volha)
        zoneStat <- zonal(volha, compsRas, "mean")
        poly[["volumeperha"]] <- zoneStat[, "mean"]
        colRamp <- colorRampPalette(c("lightgoldenrod1", "tomato2"))(10)
        polyCols <- colRamp[as.numeric(cut(poly[["volumeperha"]],
            breaks = 10))]
        plot(poly, col = polyCols, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n", main = "Volume/Ha")
        text(gCentroid(poly, byid = TRUE), round(poly[["volumeperha"]],
            1), col = "blue", font = 2)
        writeRaster(volha, "./output/volperha_mod2.tif", format = "GTiff",
            overwrite = TRUE)
        writeOGR(poly, "output", "vhamod2", driver = "ESRI Shapefile",
            overwrite = TRUE)
    }
    else if (model == "mod3") {
        mod <- lm(log(VUB) ~ I(log(DBH^2 * H)), data = data)
        coef(mod)
        vol <- function(x) {
            exp(coef(mod)[1]) * (((23.1^2) * x)^coef(mod)[2])
        }
        vol.agg <- aggregate(img, fact = 4/res(img))
        volimg <- raster::calc(vol.agg, vol)
        x <- xres(volimg)
        y <- yres(volimg)
        ha <- (x * y)/10000
        volha.func <- function(x) {
            x/ha
        }
        volha <- raster::calc(volimg, volha.func)
        par(mar = rep(0.5, 4))
        compsRas <- rasterize(poly, volha)
        zoneStat <- zonal(volha, compsRas, "mean")
        poly[["volumeperha"]] <- zoneStat[, "mean"]
        colRamp <- colorRampPalette(c("lightgoldenrod1", "tomato2"))(10)
        polyCols <- colRamp[as.numeric(cut(poly[["volumeperha"]],
            breaks = 10))]
        plot(poly, col = polyCols, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n", main = "Volume/Ha")
        text(gCentroid(poly, byid = TRUE), round(poly[["volumeperha"]],
            1), col = "blue", font = 2)
        writeRaster(volha, "./output/volperha_mod3.tif", format = "GTiff",
            overwrite = TRUE)
        writeOGR(poly, "output", "vhamod3", driver = "ESRI Shapefile",
            overwrite = TRUE)
    }
    else if (model == "mod4") {
        mod <- lm(log(VUB) ~ I(log(DBH^2)) + I(log(H^2)), data = data)
        coef(mod)
        vol <- function(x) {
            exp(coef(mod)[1]) * ((23.1^2)^coef(mod)[2]) * ((x^2)^coef(mod)[3])
        }
        vol.agg <- aggregate(img, fact = 4/res(img))
        volimg <- raster::calc(vol.agg, vol)
        x <- xres(volimg)
        y <- yres(volimg)
        ha <- (x * y)/10000
        volha.func <- function(x) {
            x/ha
        }
        volha <- raster::calc(volimg, volha.func)
        par(mar = rep(0.5, 4))
        compsRas <- rasterize(poly, volha)
        zoneStat <- zonal(volha, compsRas, "mean")
        poly[["volumeperha"]] <- zoneStat[, "mean"]
        colRamp <- colorRampPalette(c("lightgoldenrod1", "tomato2"))(10)
        polyCols <- colRamp[as.numeric(cut(poly[["volumeperha"]],
            breaks = 10))]
        plot(poly, col = polyCols, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n", main = "Volume/Ha")
        text(gCentroid(poly, byid = TRUE), round(poly[["volumeperha"]],
            1), col = "blue", font = 2)
        writeRaster(volha, "./output/volperha_mod4.tif", format = "GTiff",
            overwrite = TRUE)
        writeOGR(poly, "output", "vhamod4", driver = "ESRI Shapefile",
            overwrite = TRUE)
    }
    else if (model == "mod5") {
        mod <- lm(log(VUB) ~ I(log(H)), data = data)
        coef(mod)
        vol <- function(x) {
            exp(coef(mod)[1]) * (x^coef(mod)[2])
        }
        vol.agg <- aggregate(img, fact = 4/res(img))
        volimg <- raster::calc(vol.agg, vol)
        x <- xres(volimg)
        y <- yres(volimg)
        ha <- (x * y)/10000
        volha.func <- function(x) {
            x/ha
        }
        volha <- raster::calc(volimg, volha.func)
        par(mar = rep(0.5, 4))
        compsRas <- rasterize(poly, volha)
        zoneStat <- zonal(volha, compsRas, "mean")
        poly[["volumeperha"]] <- zoneStat[, "mean"]
        colRamp <- colorRampPalette(c("lightgoldenrod1", "tomato2"))(10)
        polyCols <- colRamp[as.numeric(cut(poly[["volumeperha"]],
            breaks = 10))]
        plot(poly, col = polyCols, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n", main = "Volume/Ha")
        text(gCentroid(poly, byid = TRUE), round(poly[["volumeperha"]],
            1), col = "blue", font = 2)
        writeRaster(volha, "./output/volperha_mod5.tif", format = "GTiff",
            overwrite = TRUE)
        writeOGR(poly, "output", "vhamod5", driver = "ESRI Shapefile",
            overwrite = TRUE)
    }
  }

wslerry/tRee documentation built on Oct. 5, 2019, 10:30 p.m.