|
29 | 29 | #' |
30 | 30 | #' @export |
31 | 31 |
|
32 | | -net_calc <- function(netdelin, vpu, nhdplus_path ){ |
| 32 | +net_calc <- function(netdelin, vpu, nhdplus_path){ |
33 | 33 | directory <- grep(paste(vpu, "/NHDPlusAttributes", sep = ""), |
34 | 34 | list.dirs(nhdplus_path, full.names = T), |
35 | 35 | value = T) |
36 | 36 | Vaa <- grep("PlusFlowlineVAA.dbf", |
37 | | - list.files(directory[1],full.names=T), |
| 37 | + list.files(directory[1], full.names = T), |
38 | 38 | value = T) |
39 | 39 | slope <- grep("elevslope.dbf", |
40 | | - list.files(directory, full.names=T), |
| 40 | + list.files(directory, full.names = T), |
41 | 41 | value = T) |
42 | | - flow.files <- grep("PlusFlow.dbf", list.files(directory[1], |
43 | | - full.names = T), value = T) |
| 42 | + flow.files <- grep("PlusFlow.dbf", |
| 43 | + list.files(directory[1], full.names = T), |
| 44 | + value = T) |
| 45 | + |
44 | 46 | flow <- foreign::read.dbf(flow.files) |
45 | 47 | vaa <- foreign::read.dbf(Vaa) |
46 | 48 | slope <- foreign::read.dbf(slope) |
47 | 49 | names(slope) <- toupper(names(slope)) |
48 | 50 | names(vaa) <- toupper(names(vaa)) |
| 51 | + |
49 | 52 | full.net <- unique(netdelin$Network) |
50 | | - reach.data <- Reduce(function(x, y) merge(x, y, |
51 | | - by.x = "net.comid", |
52 | | - by.y = "COMID", |
53 | | - all.x = T), |
54 | | - list(full.net, vaa, slope)) |
| 53 | + |
| 54 | + reach.data <- Reduce(function(x, y) |
| 55 | + merge(x, y, by.x = "net.comid", by.y = "COMID", all.x = T), |
| 56 | + list(full.net, vaa, slope)) |
| 57 | + |
55 | 58 | #calculate network order |
56 | | - WS.ord <- aggregate(reach.data[, "STREAMORDE"], |
57 | | - by = list(group.comid = reach.data[, "group.comid"]), |
58 | | - max) |
59 | | - names(WS.ord) <- c("COMID", "WS.ord") |
| 59 | + WS.ord <- reach.data[as.character(reach.data[,"group.comid"]) == |
| 60 | + as.character(reach.data[,"net.comid"]), |
| 61 | + c("net.id","M", "STREAMORDE")] |
| 62 | + |
| 63 | + names(WS.ord) <- c("net.id", "M", "WS.ord") |
| 64 | + |
60 | 65 | #catchemnts catchment area |
| 66 | + #group by, substract, multiply |
| 67 | + #value at end of flowline |
61 | 68 | cat.area <- aggregate(reach.data[, c("AREASQKM", "LENGTHKM")], |
62 | | - by = list(COMID = reach.data[, "group.comid"]), |
| 69 | + by = list(net.id = reach.data[, "net.id"], |
| 70 | + group.comid = reach.data[,"group.comid"]), |
63 | 71 | sum) |
64 | | - drain.den <- cat.area[ ,"LENGTHKM"] / cat.area[ ,"AREASQKM"] |
| 72 | + incr <- reach.data[as.character(reach.data[,"group.comid"]) == |
| 73 | + as.character(reach.data[,"net.comid"]), |
| 74 | + c("net.id", "AREASQKM","LENGTHKM", "M")] |
| 75 | + |
| 76 | + incr <- merge(incr, cat.area, by = "net.id") |
| 77 | + area <- (incr[,"AREASQKM.y"] - incr[,"AREASQKM.x"]) + incr[,"AREASQKM.x"]*incr[,"M"] |
| 78 | + len <- (incr[,"LENGTHKM.y"] - incr[,"LENGTHKM.x"]) + incr[,"LENGTHKM.x"]*incr[,"M"] |
| 79 | + |
| 80 | + #scaled length and catchment vlaues |
| 81 | + cat.area <- data.frame(net.id = incr[,"net.id"], |
| 82 | + AreaSQKM = area, LengthKM = len) |
| 83 | + |
| 84 | + drain.den <- cat.area[ ,"LengthKM"] / cat.area[ ,"AreaSQKM"] |
65 | 85 | cat.area <- data.frame(cat.area, drain.den) |
| 86 | + |
66 | 87 | #diversion feature count |
67 | 88 | #counts minor flow paths of divergences |
68 | | - if (any(reach.data[,c("STREAMORDE")] != reach.data[,"STREAMCALC"] & |
| 89 | + if (any(reach.data[,c("STREAMORDE")] != |
| 90 | + reach.data[,"STREAMCALC"] & |
69 | 91 | reach.data[,"DIVERGENCE"]==2)){ |
| 92 | + |
70 | 93 | div.rm <- reach.data[reach.data[,c("STREAMORDE")] != |
71 | | - reach.data[,"STREAMCALC"] & reach.data[, "DIVERGENCE"] == 2, |
72 | | - c("net.comid", "group.comid")] |
| 94 | + reach.data[,"STREAMCALC"] & |
| 95 | + reach.data[, "DIVERGENCE"] == 2, |
| 96 | + c("net.id", "net.comid", "group.comid")] |
73 | 97 |
|
74 | 98 | diver.cnt <- aggregate(div.rm[, "group.comid"], |
75 | | - by = list(div.rm[, "group.comid"]), |
| 99 | + by = list(div.rm[,"net.id"], div.rm[, "group.comid"]), |
76 | 100 | length) |
77 | | - names(diver.cnt) <- c("COMID", "diver.cnt") |
78 | | - } else { |
79 | | - diver.cnt<-data.frame(COMID=99999,diver.cnt=999999) |
80 | | - } |
| 101 | + |
| 102 | + names(diver.cnt) <- c("net.id", "diver.cnt") |
| 103 | + |
| 104 | + } else { |
| 105 | + diver.cnt <- data.frame(net.id = 99999, diver.cnt = 999999) |
| 106 | + } |
81 | 107 |
|
82 | 108 | #headwaters & Tribs |
83 | | - head.h2o <- aggregate(reach.data[reach.data[,"STARTFLAG"] == 1, |
84 | | - "STREAMORDE"], |
85 | | - by = list(group.comid = |
86 | | - reach.data[reach.data[,"STARTFLAG"] == |
87 | | - 1, "group.comid"]), |
88 | | - length) |
89 | | - names(head.h2o) <- c("COMID","head.h2o") |
90 | | - trib.jun <- as.numeric(as.character(head.h2o[,"head.h2o"])) - 1 |
| 109 | + head.h2o <- aggregate(reach.data[ |
| 110 | + reach.data[,"STARTFLAG"] == 1, "STREAMORDE"], |
| 111 | + by = list(reach.data[reach.data[,"STARTFLAG"] == 1, "net.id"]), |
| 112 | + length) |
| 113 | + |
| 114 | + names(head.h2o) <- c("net.id", "head.h2o") |
| 115 | + |
| 116 | + trib.jun <- as.numeric(as.character(head.h2o[, "head.h2o"])) - 1 |
91 | 117 | head.h2o <- data.frame(head.h2o, trib.jun) |
92 | 118 |
|
93 | 119 | #edge count |
94 | | - edges <- head.h2o[,"head.h2o"]+head.h2o[,"trib.jun"] |
95 | | - reach.cnt <- data.frame(COMID = head.h2o[,"COMID"], reach.cnt = edges) |
96 | | - names(reach.cnt) <- c("COMID", "reach.cnt") |
97 | | - |
98 | | - #relief |
99 | | - maxelev <- aggregate(reach.data[,"MAXELEVSMO"], |
100 | | - by = list(reach.data[,"group.comid"]), |
101 | | - max) |
102 | | - minelev <- aggregate(reach.data[, "MINELEVSMO"], |
103 | | - by = list(reach.data[, "group.comid"]), |
104 | | - min) |
105 | | - relief <- maxelev[,"x"]-minelev[,"x"] |
106 | | - relief <- data.frame(COMID = maxelev[,"Group.1"], |
107 | | - maxelev = maxelev[,"x"], |
108 | | - minelev = minelev[,"x"], |
109 | | - releif = relief) |
| 120 | + edges <- head.h2o[,"head.h2o"] + head.h2o[,"trib.jun"] |
| 121 | + reach.cnt <- data.frame(net.id = head.h2o[,"net.id"], reach.cnt = edges) |
| 122 | + |
| 123 | + #relief - at outlet; I want to move this to basin metric |
| 124 | + #maxelev <- aggregate(reach.data[,"MAXELEVSMO"], |
| 125 | + # by = list(reach.data[,"group.comid"]), |
| 126 | + # max) |
| 127 | + #minelev <- aggregate(reach.data[, "MINELEVSMO"], |
| 128 | + # by = list(reach.data[, "group.comid"]), |
| 129 | + # min) |
| 130 | + #relief <- maxelev[,"x"]-minelev[,"x"] |
| 131 | + #relief <- data.frame(COMID = maxelev[,"Group.1"], |
| 132 | + # maxelev = maxelev[,"x"], |
| 133 | + # minelev = minelev[,"x"], |
| 134 | + # releif = relief) |
110 | 135 |
|
111 | 136 | #aggregate table for summaries of group comid |
112 | | - data.out <- as.data.frame(unique(full.net[, c("group.comid", "vpu")])) |
113 | | - names(data.out)[1] <- "COMID" |
114 | | - data.out<-Reduce(function(x, y) merge(x, y, |
115 | | - by = "COMID", |
116 | | - all.x = T), |
117 | | - list(data.out, relief, head.h2o, cat.area, WS.ord, reach.cnt, diver.cnt)) |
| 137 | + data.out <- unique(full.net[, c("net.id","group.comid", "vpu")]) |
| 138 | + |
| 139 | + names(data.out)[2] <- "COMID" |
| 140 | + |
| 141 | + data.out <- Reduce(function(x, y) |
| 142 | + merge(x, y, by = "net.id", all.x = T), |
| 143 | + list(data.out, WS.ord,head.h2o, reach.cnt, diver.cnt, cat.area))#, relief)) |
| 144 | + |
118 | 145 | return(data.out) |
119 | 146 | } |
0 commit comments