# HW 3 droughtstat.R JLCY, 19 Dec 2013

############################################################ 

# find obs drought statistics

{

    # row sum of LFmonflow
    annual_flow = apply(LFmonflow, 1, sum)

    # set median as threshold
    threshold = median(annual_flow)

    max_drought = 0
    max_surplus = 0

    sum_drought = 0
    sum_surplus = 0

    # sum up total drought & surplus
    for (i in 1:nyrs) {
        diff = annual_flow[i] - threshold
        if (diff < 0) {
            sum_drought = sum_drought + diff
            if (max_drought > diff) 
                max_drought = diff
        }
        if (diff > 0) {
            sum_surplus = sum_surplus + diff
            if (max_surplus < diff) 
                max_surplus = diff
        }
    }

    cat("Max of obs drought = ", abs(max_drought), "\n")
    cat("Max of sim drought = ", max(sim_MD), ", Standard deivation = ", sd(sim_MD), 
        "\n", "\n")

    cat("Max of obs surplus = ", max_surplus, "\n")
    cat("Max of sim surplus = ", max(sim_MS), ", Standard deivation = ", sd(sim_MS), 
        "\n", "\n")

    cat("Sum of obs drought = ", abs(sum_drought), "\n")
    cat("Sum of obs surplus = ", sum_surplus, "\n", "\n")

    ############################################################ 

    where_drought = c()
    where_surplus = c()

    j = 1
    k = 1

    # locate years of drought & surplus
    for (i in 1:nyrs) {
        diff = annual_flow[i] - threshold
        if (diff < 0) {
            where_drought[j] = i
            j = j + 1
        }
        if (diff > 0) {
            where_surplus[k] = i
            k = k + 1
        }
    }

    # set first length (dummy)
    length_drought = 0
    length_surplus = 0

    j = 1
    k = 1

    # find out length of drought & surplus in years
    for (i in 2:length(where_drought)) {
        dummy = where_drought[i] - where_drought[i - 1]
        # difference > 1, new length of drought
        if (dummy > 1) {
            length_drought = c(length_drought, j)
            j = 1
        }
        # difference = 1, keep counting
        if (dummy == 1) {
            j = j + 1
        }
    }

    # similarly...
    for (i in 2:length(where_surplus)) {
        dummy = where_surplus[i] - where_surplus[i - 1]
        if (dummy > 1) {
            length_surplus = c(length_surplus, k)
            k = 1
        }
        if (dummy == 1) {
            k = k + 1
        }
    }

    # ignore the 1st length (dummy)
    length_drought = length_drought[2:length(length_drought)]
    length_surplus = length_surplus[2:length(length_surplus)]

    ############################################################ 

    # find out average, min & max length of drought & surplus

    length_drought_mean = mean(length_drought)
    length_drought_min = min(length_drought)
    length_drought_max = max(length_drought)

    cat("Average length of obs drought in years = ", length_drought_mean, "\n", 
        "\n")
    cat("Min length of obs drought in years = ", length_drought_min, "\n", "\n")
    cat("Max length of obs drought in years = ", length_drought_max, "\n")
    cat("Max length of sim drought in years = ", max(sim_LD), ", Standard deivation = ", 
        sd(sim_LD), "\n", "\n")

    length_surplus_mean = mean(length_surplus)
    length_surplus_min = min(length_surplus)
    length_surplus_max = max(length_surplus)

    cat("Average length of obs surplus in years = ", length_surplus_mean, "\n", 
        "\n")
    cat("Min length of obs surplus in years = ", length_surplus_min, "\n", "\n")
    cat("Max length of obs surplus in years = ", length_surplus_max, "\n")
    cat("Max length of sim surplus in years = ", max(sim_LS), ", Standard deivation = ", 
        sd(sim_LS), "\n", "\n")

}

## Error: 找不到物件 'LFmonflow'