Set 3 PrepareWetlandDataForStreamCatAccumulation.Rmd

---
title: "Prep Wetland Data for StreamCat Accumulation"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  html_document:
    theme: yeti
    highlighted: default 
    toc: yes
    toc_float: true
---


## Code to prepare data for StreamCat processing

### Calculate wetland area by **wetlands type** for StreamCat accumulation

```{r, eval=F}
library(foreign)
library(dplyr)

year='2011'

working_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/'
rds_path <- paste0(working_dir, '/Wetlands_NLCD',year,'/FinalTables/')
table_path = paste0(rds_path, 'WetlandTables/')
zonal_path <- paste0(working_dir, 'Wetlands_NLCD2011/WetlandCat/Zonal/')
nhd_path <- 'D:/GISData/NHDPlusV21'
n_zonal_path <- paste0(working_dir, 'Wetlands_NLCD2011/NitrogenModeling/Allocation_Accumulation/')

fullwetlands <- readRDS(paste0(rds_path, 'WetConnectMetrics_',year,'.rds'))
fullwetlands$travel <- ifelse(fullwetlands$Type == 'NRSur', fullwetlands$MagOv, fullwetlands$MagSh)
fullwetlands$travel[fullwetlands$Type == 'Ripar'] <- 0
fullwetlands$travel <- fullwetlands$travel + 1
fullwetlands <- subset(fullwetlands, select = c(WetId,COMID,Type2,travel))
fullwetlands$wetid_comid <- paste0(fullwetlands$WetId, '_', fullwetlands$COMID)

rpucomid <- read.csv('L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/COMID_HydroRegion_RPU.csv')[c('COMID','RPU')]

hydroregions <- read.csv(paste0(working_dir, 'hydro-regions.csv'))
names(hydroregions) <- c('Region', 'RPU', 'Hydro')

hydro.rgns <- c("01","02","03S","03N","03W","04","05","06","07","08","09","10L","10U","11","12","13","14","15","16","17","18")
rgns <- c("NE","MA","SA","SA","SA","GL","MS","MS","MS","MS","SR","MS","MS","MS","TX","RG","CO","CO","GB","PN","CA")

for(i in 1:length(hydro.rgns)){
  print(hydro.rgns[i])
  
  cats = read.dbf(paste0(nhd_path,"/NHDPlus",rgns[i],"/NHDPlus",hydro.rgns[i],"/NHDPlusCatchment/Catchment.dbf"))[, c('GRIDCODE','FEATUREID')]
  
  rpus <- hydroregions$RPU[hydroregions$Hydro %in% hydro.rgns[i]]
  
  #Combine wetland basin area file from RPUs for each VPU
  #Any of the final tables could be used. Simply need COMID-WETLAND ID connection
  for(j in 1:length(rpus)){
    tmp <- read.csv(paste0(table_path, 'BasinAreaSqKm_', rpus[j], '.csv'))[, c('WET_ID','COMID','CatAreaSqKm','WetlandAreaSqKm')]
    tmp <- merge(tmp, rpucomid, by='COMID')
    tmp <- tmp[tmp$RPU %in% rpus[j], ]
    if(j == 1){
      zonal <- tmp
    }else{
      zonal <- rbind(zonal, tmp)
    }
  }
  
  zonal <- zonal[!duplicated(zonal), ]
  zonal$wetid_comid <- paste0(zonal$WET_ID, '_', zonal$COMID)
  zonal <- subset(zonal, select=-COMID)
  
  zonal <- merge(zonal, fullwetlands, by='wetid_comid')

  zonal$SUM <- zonal$WetlandAreaSqKm * 100
  zonal$COUNT <- (zonal$WetlandAreaSqKm * 1e6) / 900

  zonal0 <- zonal %>% group_by(Type2, COMID) %>% summarize(SUM=sum(SUM),
                                                          COUNT=sum(COUNT),
                                                          AREA=max(CatAreaSqKm)*1e6)
  
  zonal0 <- merge(zonal0, cats, by.x='COMID', by.y='FEATUREID')
  names(zonal0)[length(zonal0)] <- 'VALUE'
  zonal0 <- zonal0[, c('VALUE','COUNT','AREA','SUM','Type2')]
  
  wet_types <- unique(na.omit(zonal$Type2))
  
  for(j in 1:length(wet_types)){
    z0 <- na.omit(zonal0[zonal0$Type2 == wet_types[j], ])
    write.dbf(z0, file=paste0(n_zonal_path, 'DBF_stash/', 'WtAreasByType2_', wet_types[j], '_', hydro.rgns[i], '.dbf'))
  }
}
```


## Aggregating StreamCat output into long files

### Combine accumulated wetland areas by wetland type

```{r, eval=F}
accum_dir <- "L:/Priv/CORFiles/Geospatial_Library_Projects/WetlandConnectivity/SpatialDataInputs/Wetlands_NLCD2011b/NitrogenModeling/Allocation_Accumulation/"
nmodel_dir <- paste0(working_dir,'nitrogen-modeling/')
Types <- c("Ripar","NRMid","NRDeep","NRShallow")
hydro.rgns <- c("01","02","03S","03N","03W","04","05","06","07","08","09","10L","10U","11","12","13","14","15","16","17","18")

for(i in 1:length(hydro.rgns)){
  print(hydro.rgns[i])
  for(k in 1:length(Types)){
    print(Types[k]) 
    tmp <- read.csv(paste0(accum_dir, 'WtAreasByType2_', Types[k], '_', hydro.rgns[i], '.csv'))
    #names(tmp)[2] <- paste0(Types[k], 'WtUpCount')
    if(k == 1){
      tmp <- tmp[, c('COMID','UpCatAreaSqKm','UpCatCount')]
      names(tmp)[3] <- paste0(Types[k], 'WtUpCount')
      tmpdf <- tmp
    }else{
      tmp <- tmp[, c('COMID','UpCatCount')]
      names(tmp)[2] <- paste0(Types[k], 'WtUpCount')
      tmpdf <- merge(tmpdf, tmp, by='COMID')
    }
  }
  if(i == 1){
    outdf <- tmpdf
  }else{
    outdf <- rbind(outdf, tmpdf)
  }
}

saveRDS(outdf, paste0(nmodel_dir, 'splitcat-metrics/Wetland2PxlUpCount.rds'))
#saveRDS(outdf, paste0(nmodel_dir, 'splitcat-metrics/WetlandPxlUpCount.rds'))

```


## Prepare SplitCat metrics for combining with upstream StreamCat metrics

### Calculate N applied to NRSA SplitCats (entire area)

```{r, eval=F}
library(raster)

sc_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/NARS_SplitCat_Metrics/NRSA-2008-09/split_cat/splits/split_full/'
translation_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/NARS_SplitCat_Metrics/NRSA-2008-09/'
nraster_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/WetlandConnectivity/NitrogenModeling/LandscapeRasters/'
cmaq_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/LandscapeRasters/QAComplete/'
wet_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/Wetlands_NLCD2011/AllWetlands/'
wetcat_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/Wetlands_NLCD2011/WetlandCat/WetCats/'
nmodel_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/nitrogen-modeling/'
accum_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/Allocation_and_Accumulation/'

ndata <- read.csv(paste0(nmodel_dir, 'N.data.csv'))
translation <- read.csv(paste0(translation_dir, 'sites.csv'))[, c('SITE_ID','UNIQUE_ID')]
ndata <- merge(ndata, translation, by='SITE_ID')
cbnf <-  raster(paste0(nraster_dir, 'cbnf.tif'))
manure <- raster(paste0(nraster_dir, 'manure.tif'))
fert <- raster(paste0(nraster_dir, 'fert.tif'))
cmaq <- raster(paste0(cmaq_dir, 'TD_N.tif'))

tokg <- function(inras){
  count <- sum(!is.na(inras@data@values))
  sum <- cellStats(inras, 'sum')
  ha <- count * 900 * 0.0001 
  kg <- (sum / count) * ha
  return(kg)
}

sites <- as.character(ndata$UNIQUE_ID)
siteids <- as.character(ndata$SITE_ID)
for(i in 1:length(sites)){
  print(sites[i])
  sc <- paste0(sc_dir, sites[i], '.tif')
  if(file.exists(sc)){
    sc <- raster(sc)
    sc[!is.na(sc)] <- 1
    
    n1 <- crop(cbnf, extent(sc)) * sc
    n2 <- crop(manure, extent(sc)) * sc
    n3 <- crop(fert, extent(sc)) * sc
    n4 <- crop(cmaq, extent(sc)+c(-10000,10000,-10000,10000))
    n4 <- resample(n4, sc, method="ngb") * sc
    SUM <- tokg(n1) + tokg(n2) + tokg (n3) + tokg(n4)
    COUNT <- sum(!is.na(sc@data@values))
    if(i == 1){
      outdf <- data.frame(SITE_ID = siteids[i], UNIQUE_ID=sites[i], SUM=SUM, COUNT=COUNT)  
    }else{
      outdf <- rbind(outdf, 
                     data.frame(SITE_ID = siteids[i], UNIQUE_ID=sites[i], SUM=SUM, COUNT=COUNT))
    }
  }
}

tmp <- ndata[, c('SITE_ID','COMID')]
outdf <- merge(outdf, tmp, by='SITE_ID')
write.csv(outdf, paste0(nmodel_dir, 'splitcat-fullcatchment/full-splitcat-metrics.csv'), row.names = F)

#tester <- read.csv(paste0(nmodel_dir, 'splitcat-fullcatchment/full-splitcat-metrics.csv'))


```


### Calculate N applied to **each wetland catchment** (by IDs) within NRSA SplitCats

TABLE INFO:
SITE_ID = NRSA SPLITCAT ID, WET_ID = WETLAND ID, WtCatCount = NUM. WETLAND CATCHMENT PIXELS
SUM = SUM TOTAL N (KG), CatCount = NUM. SPLITCAT PIXELS
rpu = NHDPLUS RASTER PROCESSING UNIT, COMID = NHDPLUS UNIQUE ID

* It is possible to have wetland areas that are larger than wetland catchment areas because of where NAs from full streams land and where the edge of splitcats occur. These have been verified to be correct.


```{r, eval=F}
library(raster)
library(dplyr)
# Define folders
sc_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/split_full/'
translation_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/NARS_SplitCat_Metrics/NRSA-2008-09/'
#nraster_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/WetlandConnectivity/NitrogenModeling/LandscapeRasters/'
#cmaq_dir <- 'L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/LandscapeRasters/QAComplete/'
nmodel_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/nitrogen-modeling/'
wetcat_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/Wetlands_NLCD2011/WetlandCat/WetCats/'
wetland_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/Wetlands_NLCD2011/AllWetlands_rpu/'
rpucomid <- read.csv('L:/Priv/CORFiles/Geospatial_Library_Projects/StreamCat/COMID_HydroRegion_RPU.csv')[c('COMID','RPU')]
ndata <- read.csv(paste0(nmodel_dir, 'N.data.csv'))
ndata <- merge(ndata, rpucomid, by='COMID')
translation <- read.csv(paste0(translation_dir, 'sites.csv'))[, c('SITE_ID','UNIQUE_ID')]
ndata <- merge(ndata, translation, by='SITE_ID')
# Read in rasters
cbnf <-  raster(paste0(nraster_dir, 'cbnf.tif'))
manure <- raster(paste0(nraster_dir, 'manure.tif'))
fert <- raster(paste0(nraster_dir, 'fert.tif'))
cmaq <- raster(paste0(cmaq_dir, 'TD_N.tif'))
#Create zonal stats function to get N to each wetland catchment in each splitcat
zonalstat <- function(inras, zonalras){
  inzone <- zonalras@data@values
  yes <- !is.na(inzone)
  if(sum(yes) != 0){
    ras <- inras@data@values[yes]
    inzone <- inzone[yes]
    tmp <- data.frame(ras=ras, inzone=inzone)
    tmp <- data.frame(tmp %>% group_by(inzone) %>% summarise(SUM=sum(ras), COUNT=n()))
    tmp$ha <- tmp$COUNT * 900 * 0.0001
    tmp$kg <- (tmp$SUM / tmp$COUNT) * tmp$ha
    tmp$SUM <- tmp$kg
    names(tmp)[1] <- 'WET_ID'
    #tmp$COUNT <- sum(yes)
    tmp <- tmp[ , c('WET_ID','COUNT','SUM')]
  }else{
    tmp <- data.frame(WET_ID=NA, SUM=0, COUNT=0)
  }
  return(tmp)
}

zonecount <- function(zonalras){
  inzone <- zonalras@data@values
  yes <- !is.na(inzone)
  if(sum(yes) != 0){
    inzone <- inzone[yes]
    tmp <- data.frame(inzone=inzone)
    tmp <- data.frame(tmp %>% group_by(inzone) %>% summarise(COUNT=n()))
    names(tmp)[1] <- 'WET_ID'
  }else{
    tmp <- data.frame(WET_ID=NA, COUNT=0)
  }
  return(tmp)
}

#Loop through RPUs to get wetland catchments
allsites <- as.character(ndata$SITE_ID)
unqids <- as.character(ndata$UNIQUE_ID)
rpus <- as.character(ndata$RPU)
unqrpus <- unique(rpus)
unqrpus <- unqrpus[order(unqrpus)]
for(k in 1:length(unqrpus)){
  print(unqrpus[k])
  wetcat <- raster(paste0(wetcat_dir, 'WetlandCat_', unqrpus[k], '.tif'))
  wetarea <- raster(paste0(wetland_dir, 'WetlandsRgnGrp_', unqrpus[k], '.tif'))
  sites <- allsites[rpus == unqrpus[k]]
  uniques <- unqids[rpus == unqrpus[k]]
  #Loop through splitcats and calc total kg of N
  for(i in 1:length(uniques)){
    print(uniques[i])
    sc <- paste0(sc_dir, uniques[i], '.tif')
    if(file.exists(sc)){
      sc <- raster(sc)
      #sc[!is.na(sc)] <- 1
      #read and crop rasters by splitcat
      #wc <- crop(wetcat, extent(sc)) * sc
      if(!all(is.na(values(wc)))){
        wt <- crop(wetarea, extent(sc)) * sc
        # n1 <- crop(cbnf, extent(sc)) * sc
        # n2 <- crop(manure, extent(sc)) * sc
        # n3 <- crop(fert, extent(sc)) * sc
        # n4 <- crop(cmaq, extent(sc)+c(-10000,10000,-10000,10000))
        # n4 <- resample(n4, sc, method="ngb") * sc
        #Cal zonal values for each wetland
        #tmp <- zonalstat(n1,wc)
        #tmp$SUM <- tmp$SUM + zonalstat(n2,wc)$SUM + zonalstat(n3,wc)$SUM + zonalstat(n4,wc)$SUM
        tmp <- zonecount(wt)
        names(tmp)[2] <- 'WtCount'
        tmp$SITE_ID <- sites[i]
        tmp$UNIQUE_ID <- uniques[i]
        tmp$CatCount <- cellStats(sc, sum)
        tmp$rpu <- unqrpus[k]
      }else{
        tmp <- data.frame(WET_ID=NA,  WtCount=0, 
                          SITE_ID=sites[i], UNIQUE_ID=uniques[i], 
                          CatCount=cellStats(sc, sum), 
                          rpu=unqrpus[k])
      }
    }else{
      tmp <- data.frame(WET_ID=NA, WtCount=NA, 
                        SITE_ID=sites[i], UNIQUE_ID=uniques[i], 
                        CatCount=NA, 
                        rpu=unqrpus[k])
    }
    if(i == 1){
      tmpdf <- tmp
    }else{
      tmpdf <- rbind(tmpdf, tmp)
    }
  }
  if(k == 1){
    outdf <- tmpdf
  }else{
    outdf <- rbind(outdf, tmpdf)
  }
}

# Pair with NHDPlus COMID
write.csv(outdf, paste0(nmodel_dir, 'splitcat-wtcount-metrics.csv'), row.names = F)
```

### Combine splitcat wtcount with upcatchment wtcounts 

```{r, eval=F}
library(dplyr);library(tidyr)
rds_path <- paste0(working_dir, '/Wetlands_NLCD',year,'/FinalTables/')
nmodel_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/nitrogen-modeling/'
ndata <- read.csv(paste0(nmodel_dir, 'N.data.csv'))
ndata <- merge(ndata, rpucomid, by='COMID')
translation <- read.csv(paste0(translation_dir, 'sites.csv'))[, c('SITE_ID','UNIQUE_ID')]
ndata <- merge(ndata, translation, by='SITE_ID')
ndata <- ndata[, c('SITE_ID','UNIQUE_ID','COMID')]

fullwetlands <- readRDS(paste0(rds_path, 'WetConnectMetrics_',year,'.rds'))
fullwetlands <- fullwetlands[, c('WetId','Type2')]
names(fullwetlands) <- c('WET_ID','Type2')
spc <- read.csv(paste0(nmodel_dir, 'splitcat-wtcount-metrics.csv'))
spc <- merge(spc, fullwetlands, by='WET_ID', all.x=T)
nowet <- spc[is.na(spc$WET_ID), ]

spc <- spc %>% group_by(SITE_ID, Type2) %>% summarise(WtCount=sum(WtCount),
                                             CatCount=max(CatCount))
spc <- spread(spc, Type2, WtCount)
spc <- spc[, -ncol(spc)]

ndata <- merge(ndata, spc, by='SITE_ID', all.x=T)
no_catcount <- ndata[is.na(ndata$CatCount), ]
ndata <- ndata[!is.na(ndata$CatCount), ]


upcounts <- readRDS(paste0(nmodel_dir, 'splitcat-metrics/Wetland2PxlUpCount.rds'))

ndata <- merge(ndata, upcounts, by='COMID')
ndata[is.na(ndata)] <- 0

ndata$WsArea <- ((ndata$CatCount * 900) / 1e6) + ndata$UpCatAreaSqKm
ndata$WsNRDeep <- ((ndata$NRDeep * 900) / 1e6) + ((ndata$NRDeepWtUpCount*900) / 1e6)
ndata$WsNRMid <- ((ndata$NRMid * 900) / 1e6) + ((ndata$NRMidWtUpCount*900) / 1e6)
ndata$WsNRShallow <- ((ndata$NRShallow * 900) / 1e6) + ((ndata$NRShallowWtUpCount*900) / 1e6)
ndata$WsRipar <- ((ndata$Ripar * 900) / 1e6) + ((ndata$RiparWtUpCount*900) / 1e6)

ndata <- ndata[, c('COMID','SITE_ID','UNIQUE_ID','WsArea','WsNRDeep','WsNRMid','WsNRShallow','WsRipar')]
ndata$pctNRDeep <- (ndata$WsNRDeep / ndata$WsArea) * 100
ndata$pctNRMid <- (ndata$WsNRMid / ndata$WsArea) * 100
ndata$pctNRShallow <- (ndata$WsNRShallow / ndata$WsArea) * 100
ndata$pctRipar <- (ndata$WsRipar / ndata$WsArea) * 100

saveRDS(ndata, paste0(nmodel_dir,'chem-modleing-new-classes-2021.04.25.rds'))

```


### Summarize wetland catchment N for each SITE_ID and by wetland type for each SITE_ID

```{r, eval=F}
library(dplyr)
year='2011'
rds_path <- paste0(working_dir, '/Wetlands_NLCD',year,'/FinalTables/')
nmodel_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/nitrogen-modeling/'

fullwetlands <- readRDS(paste0(rds_path, 'WetConnectMetrics_',year,'.rds'))

fullwetlands$travel <- ifelse(fullwetlands$Type == 'NRSur', fullwetlands$MagOv, fullwetlands$MagSh)
fullwetlands$travel[fullwetlands$Type == 'Ripar'] <- 0
fullwetlands$travel <- fullwetlands$travel + 1
fullwetlands <- subset(fullwetlands, select = c(WetId,COMID,Type2,travel,WetAreaSqKm))
fullwetlands$wetidcomid <- paste0(fullwetlands$WetId, '_', fullwetlands$COMID)

zonal <- read.csv(paste0(nmodel_dir, 'splitcat-wetlandbasin/splitcat-wetland-basin-metrics.csv'))
NAs <- zonal[is.na(zonal$WET_ID), ]
NAs <- subset(NAs, select=c(SITE_ID,SUM,WtCatCount,CatCount,WtCount,COMID))
names(NAs)[2] <- 'CatSum'
write.csv(NAs, paste0(nmodel_dir, 'splitcat-wetlandbasin/splitcat-wetland-basin-metrics--NAs.csv'), row.names=F)
zonal <- zonal[!is.na(zonal$WET_ID), ]
zonal$wetidcomid <- paste0(zonal$WET_ID, '_', zonal$COMID)
zonal <- subset(zonal, select = -COMID)
zonal <- merge(zonal, fullwetlands, by='wetidcomid')
zonal <- subset(zonal, select = -wetidcomid)
zonal$SUM[is.na(zonal$SUM)] <- 0
#SUM already as Kg
zonal$kg <- zonal$SUM
zonal$CatAreaSqKm <- (zonal$CatCount * 900) / 1e6
zonal$WetAreaSqKm <- (zonal$WtCount * 900) / 1e6
zonal$newsum <- zonal$kg / (zonal$WetAreaSqKm * 100)
zonal$newsum2 <- zonal$kg / (zonal$travel)
zonal$newsum3 <- zonal$kg / (zonal$WetAreaSqKm * 100) / (zonal$travel)

zonal0 <- zonal %>% group_by(SITE_ID) %>% summarize(CatSum=sum(kg),
                                                    WtCatCount=sum(WtCatCount),
                                                    CatCount=max(CatCount),
                                                    WtCount=sum(WtCount),
                                                    COMID=max(COMID))
zonal0 <- rbind(zonal0, NAs)

zonal1 <- zonal %>% group_by(SITE_ID) %>% summarize(CatSum=sum(newsum),
                                                    WtCatCount=sum(WtCatCount),
                                                    CatCount=max(CatCount),
                                                    WtCount=sum(WtCount),
                                                    COMID=max(COMID))
zonal1 <- rbind(zonal1, NAs)

zonal2 <- zonal %>% group_by(SITE_ID) %>% summarize(CatSum=sum(newsum2),
                                                    WtCatCount=sum(WtCatCount),
                                                    CatCount=max(CatCount),
                                                    WtCount=sum(WtCount),
                                                    COMID=max(COMID))
zonal2 <- rbind(zonal2, NAs)

zonal3 <- zonal %>% group_by(SITE_ID) %>% summarize(CatSum=sum(newsum3),
                                                    WtCatCount=sum(WtCatCount),
                                                    CatCount=max(CatCount),
                                                    WtCount=sum(WtCount),
                                                    COMID=max(COMID))
zonal3 <- rbind(zonal3, NAs)

write.csv(zonal0, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_TotalKg_splitcat.csv'), row.names = F)
write.csv(zonal1, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_ByWtArea_splitcat.csv'), row.names = F)
write.csv(zonal2, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_ByTravelTime_splitcat.csv'), row.names = F)
write.csv(zonal3, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_ByAreaTime_splitcat.csv'), row.names = F)


#Now summarize by wetland type
NAs$Type2 <- NA
NAs <- NAs[, c('Type2','SITE_ID','CatSum','WtCatCount','CatCount','WtCount','COMID')]
NAs$Type2 <- as.character(NAs$Type2)


zonal0 <- zonal %>% group_by(Type2, SITE_ID) %>% summarize(CatSum=sum(kg),
                                                          WtCatCount=sum(WtCatCount),
                                                          CatCount=max(CatCount),
                                                          WtCount=sum(WtCount),
                                                          COMID=max(COMID))
zonal0 <- rbind(as.data.frame(zonal0), NAs)

zonal1 <- zonal %>% group_by(Type2, SITE_ID) %>% summarize(CatSum=sum(newsum),
                                                          WtCatCount=sum(WtCatCount),
                                                          CatCount=max(CatCount),
                                                          WtCount=sum(WtCount),
                                                          COMID=max(COMID))
zonal1 <- rbind(as.data.frame(zonal1), NAs)

zonal2 <- zonal %>% group_by(Type2, SITE_ID) %>% summarize(CatSum=sum(newsum2),
                                                          WtCatCount=sum(WtCatCount),
                                                          CatCount=max(CatCount),
                                                          WtCount=sum(WtCount),
                                                          COMID=max(COMID))
zonal2 <- rbind(as.data.frame(zonal2), NAs)

zonal3 <- zonal %>% group_by(Type2, SITE_ID) %>% summarize(CatSum=sum(newsum3),
                                                          WtCatCount=sum(WtCatCount),
                                                          CatCount=max(CatCount),
                                                          WtCount=sum(WtCount),
                                                          COMID=max(COMID))
zonal3 <- rbind(as.data.frame(zonal3), NAs)

write.csv(zonal0, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_TotalKg_ByType_splitcat.csv'), row.names = F)
write.csv(zonal1, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_ByWtArea_ByType_splitcat.csv'), row.names = F)
write.csv(zonal2, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_ByTravelTime_ByType_splitcat.csv'), row.names = F)
write.csv(zonal3, file=paste0(nmodel_dir, 'splitcat-metrics/', 'TN_ByAreaTime_ByType_splitcat.csv'), row.names = F)

```

## Combine wetland basin splitcat N (BY TYPE) with StreamCat wetland accumulations to create site-specific summaries

```{r, eval=F}
accum_dir <- "L:/Priv/CORFiles/Geospatial_Library_Projects/WetlandConnectivity/SpatialDataInputs/Wetlands_NLCD2011b/NitrogenModeling/Allocation_Accumulation/"
working_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/'
nmodel_dir <- paste0(working_dir,'nitrogen-modeling/')

hydro.rgns <- c("01","02","03S","03N","03W","04","05","06","07","08","09","10L","10U","11","12","13","14","15","16","17","18")

#sc <- read.csv(paste0(nmodel_dir, 'splitcat-metrics/', 'TN_TotalKg_splitcat.csv'))
Types <- c("Ripar","NRMid","NRDeep","NRShallow")

for(k in 1:length(hydro.rgns)){
  print(hydro.rgns[k])
  for(i in 1:length(Types)){
    print(Types[i])
    tmp <- read.csv(paste0(accum_dir, 'TN_TotalKg_', Types[i],  '_', hydro.rgns[k], '.csv'))
    if(i == 1){
      tmp <- tmp[, c('COMID', 'UpCatAreaSqKm', 'UpCatCount', 'UpCatSum')]
      #tmp$UpCatSum <- tmp$UpCatSum / (900/10000)
      names(tmp) <- c('COMID', 'UpCatAreaSqKm', paste0(Types[i], c('UpWtCatCount','UpWtCatSum')))
      out1 <- tmp
    }else{
      tmp <- tmp[, c('COMID', 'UpCatCount', 'UpCatSum')]
      #tmp$UpCatSum <- tmp$UpCatSum / (900/10000)
      names(tmp) <- c('COMID', paste0(Types[i], c('UpWtCatCount','UpWtCatSum')))
      out1 <- merge(out1, tmp, by='COMID')
    }
  }
  
  if(k == 1){
    outdf <- out1
  }else{
    outdf <- rbind(outdf, out1)
  }
}

outdf <- outdf[!duplicated(outdf), ]
saveRDS(outdf, paste0(nmodel_dir, 'splitcat-metrics/TN_TotalKg_UpCatMetrics_ByWetType.rds'))
```

## Redo of wetland types 2021.04.24

```{r, eval=F}
working_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/'
nmodel_dir <- paste0(working_dir,'nitrogen-modeling/')



wtpxupcount <- readRDS(paste0(nmodel_dir, 'splitcat-metrics/Wetland2PxlUpCount.rds'))


```

## Combine full and wetland-filtered estimates of watershed N 

```{r, eval=F}
working_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/'
nmodel_dir <- paste0(working_dir,'nitrogen-modeling/')

upwtcatpxbytype <- readRDS(paste0(nmodel_dir, 'splitcat-metrics/WetCatPxlUpCount.rds'))
upwtcat <- readRDS(paste0(nmodel_dir, 'splitcat-metrics/TN_TotalKg_UpCatMetrics_ByWetType.rds'))
upwtcat <- upwtcat[, c('COMID', grep('UpWtCatSum', names(upwtcat), value=T))]
names(upwtcat) <- c('COMID', 'RiparUpSum', 'NRSurUpSum', 'NRSubPdUpSum', 'NRSubWdUpSum')
dwnslp <- readRDS(paste0(nmodel_dir, 'splitcat-metrics/CountDwnSlpWtWsMetrics.rds'))
wtpxupcount <- readRDS(paste0(nmodel_dir, 'splitcat-metrics/Wetland2PxlUpCount.rds'))

sc <- read.csv(paste0(nmodel_dir, 'splitcat-metrics/TN_TotalKg_ByType_splitcat.csv'))
scfull <- read.csv(paste0(nmodel_dir, 'splitcat-metrics/TN_TotalKg_splitcat.csv'))[, c('SITE_ID','CatCount')]

nout <- read.csv(paste0(nmodel_dir, 'N.data.csv'))
#nout <- nout[!(nout$SITE_ID %in% c('FW08AR218','FW08WI083')),] #Problem splitcat delineations

nout <- subset(nout, select=c(SITE_ID, COMID))

Types <- c('Ripar','NRSur','NRSubPd','NRSubWd')
#Take long tables of wetland type by COMID
for(i in 1:length(Types)){
  print(Types[i])
  tmp <- sc[sc$Type == Types[i], ]
  tmp <- tmp[, c('SITE_ID','CatSum','WtCatCount','WtCount')]
  nms <- paste0(Types[i], c('CatSum','WtCatCount','WtCount'))
  names(tmp) <- c('SITE_ID', nms)
  nout <- merge(nout, tmp, by='SITE_ID', all.x=T)
}

nout <- merge(nout, wtpxupcount, by='COMID')
nout <- merge(nout, upwtcatpxbytype, by='COMID')
nout <- merge(nout, upwtcat, by='COMID')
nout[is.na(nout)] <- 0
nout <- merge(nout, scfull, by='SITE_ID')
nout <- merge(nout, dwnslp, by='COMID')
nout$DSlpWtWs[is.na(nout$DSlpWtWs)] <- -1
nout$DSlpWtNrWs[is.na(nout$DSlpWtNrWs)] <- -1

write.csv(nout, paste0(nmodel_dir, 'splitcat-metrics/fullwatershed-splitcat-wetland-metrics-BYWETLAND-TYPE.csv'), row.names = F)

```

## Final data prep for modeling

```{r, eval=F}
nmodel_dir <- 'D:/WorkFolder/WetlandConnectivity/WetConnect_Oct2019/nitrogen-modeling/'

ninput <- read.csv(paste0(nmodel_dir, 'splitcat-metrics/fullwatershed-splitcat-wetland-metrics-BYWETLAND-TYPE.csv'))
nout <- read.csv(paste0(nmodel_dir, 'N.data.csv'))
fullWsN <- read.csv(paste0(nmodel_dir, 'splitcat-fullcatchment/fullwatershed-splitcat-metrics.csv')) [, c('SITE_ID','KgNWs')]
names(fullWsN)[2] <- 'AllKgNWs'
ninput <- merge(ninput, fullWsN, by='SITE_ID')

#Prep N export data
#nout <- nout[, c('SITE_ID','TN','TP','DOC','q.ann')]
nout <- nout[, c('SITE_ID','TN','q.ann')]
#TN(mgN/L)*q.ann(m^3/s)*sec*min*hour*days*L/m^3 / mg/Kg = Kg/yr
nout$TNKgYr <- (nout$TN * nout$q.ann * 60 * 60 * 24 * 365.24 * 1000) / 1e6

#Prep N input data
sums <- grep('Sum', names(ninput), value = T)
ninput$WtKgNWs <- rowSums(ninput[sums])
#counts <- grep('WtCatCount', names(ninput), value=T)
ninput$WsAreaHa <- (((ninput$CatCount * 900) / 1e6) * 100) + (ninput$UpCatAreaSqKm * 100)

Types <- c('Ripar','NRSur','NRSubPd','NRSubWd')
for(i in 1:length(Types)){
  wt <- paste0(Types[i], c('WtCount','WtUpCount'))
  #print(head(ninput[,wt]))
  ninput$tmpcol <- rowSums(ninput[,wt])
  ninput$tmpcol <- (ninput$tmpcol * 900) / 10000
  names(ninput)[length(ninput)] <- paste0(Types[i], 'WtAreaHa')
  
  wt <- paste0(Types[i], c('CatSum','UpSum'))
  #print(head(ninput[,wt]))
  ninput$tmpcol <- rowSums(ninput[,wt])
  names(ninput)[length(ninput)] <- paste0(Types[i], 'WtKgN')
  
  wt <- paste0(Types[i], c('WtCatCount','UpWtCatCount'))
  #print(head(ninput[,wt]))
  ninput$tmpcol <- rowSums(ninput[,wt])
  ninput$tmpcol <- (ninput$tmpcol * 900) / 10000
  names(ninput)[length(ninput)] <- paste0(Types[i], 'WtCatAreaHa')
}


ninput$WtCatAreaHaWs <- (rowSums(ninput[grep('WtCatCount', names(ninput))]) * 900) / 10000
ninput <- ninput[, c('SITE_ID','COMID', 'DSlpWtWs', 'DSlpWtNrWs',
                     grep('KgN', names(ninput), value=T),
                     grep('AreaHa', names(ninput), value=T))]

ndata <- merge(nout, ninput, by='SITE_ID')
ndata$NwKgNWs <- ndata$AllKgNWs - ndata$WtKgNWs
#ndata$q.ann_ha <- ndata$q.ann / ndata$WsAreaHa

write.csv(ndata, paste0(nmodel_dir, 'splitcat-metrics/nmodel-data.csv'), row.names = F)

```