AAv15 analysis

Restarting with the addition of ILNH data

Introduction

Objectives

1. Compare the number and area of wetlands mapped by 3-factor and NWI wetlands across study sites

2. Evaluate correlates of differences in mapping outcomes

Methods

Map comparison metrics

Total area mapped:

The simple sum of all wetland polygons for a given mapping protocol.

For 3-factor wetlands: \[{A_{3f}} = \sum_{i=1}^{n}A_{i}\]

For NWI wetlands: \[{A_{NWI}} = \sum_{i=1}^{n}A_{i}\]

Total polygons mapped:

The count of all wetland polygons for a given mapping protocol.

The formulation for 3-factor wetlands: \[{n_{3f poly}} = \sum_{i=1}^{n}n_{i}\]

The formulation for NWI wetlands: \[{A_{NWI poly}} = \sum_{i=1}^{n}n_{i}\]

Normalized Difference Area Index:

\[NDAI = \frac{({A_{3f}} - {A_{NWI}})}{({A_{3f}} + A_{NWI})}\]

Goodness of Fit [add citation]

\[G_{fit} = \frac{({A_{union}})}{({A_{3f.only}} + A_{union})} * \frac{({A_{union}})}{({A_{NWI.only}} + A_{union})}\]

Results

Assessment Area characterstics

Total area and number of AA

OGR data source with driver: ESRI Shapefile 
Source: "data/Final_shapefile/v15_data", layer: "polyAAidv15fin_centroid"
with 1751 features
It has 8 fields
[1] "OBJECTID"   "AAv14id"    "Shape_Leng" "Shape_Area" "AAidAcre"  
[6] "AAareaHa"   "CenLat"     "CenLon"    

## add the lookup table...
v15_centroid_EcoregDistDiv_LU <- read_csv("D:/Dropbox/PROJECTS/CORPS_NWI3P/RStudio_Projects/NWI_Analysis2016/data/Final_shapefile/V15_data/v15_centroid_EcoregDistDiv_LU.csv")
# and add the nwi meta lookup
## import the NWI metadata
AAv15centrNWImeta_LU <- read_csv("D:/Dropbox/PROJECTS/CORPS_NWI3P/RStudio_Projects/NWI_Analysis2016/data/Final_shapefile/V15_data/AAv15centrNWImeta_LU.csv")
nwi.meta <- AAv15centrNWImeta_LU %>% 
  # names()
  dplyr::select(c(AAv14id,IMAGE_YR, EMULSION))
#---------
# View(v15_centroid_EcoregDistDiv_LU)
# v15_centroid_EcoregDistDiv_LU %>%
#   names() %>%
#   as_tibble()
# unique(v15_centroid_EcoregDistDiv_LU$STATE)
## select the fields I want
ecoreg.lu <- v15_centroid_EcoregDistDiv_LU %>% 
  dplyr::select(AAv14id, 
District,
DIVISION,
US_L4CODE,              
US_L4NAME,          
US_L3NAME,              
NA_L3NAME,              
NA_L2NAME,
BlyDOMAIN,              
BlyDIVISION,
PROVINCE,               
SECTION,
STATE,
isoBioCli
)
  
# datatable(ecoreg.lu)
## Left join the ecoregion etc. lu with the AA
aaCtrd.ecoreg <- left_join(aacentr.tbl, ecoreg.lu, by = "AAv14id")
# Now join the NWI metadata 
aaCtrd.ecoreg <- left_join(aaCtrd.ecoreg, nwi.meta, by = "AAv14id")
## just a check to see I wave one poly per AAv14id
# aaCtrd.ecoreg %>% 
#   group_by(AAv14id) %>%
#   summarize(cnt= n()) %>% 
#   filter(cnt==1)

Summary datatabe for AA

aaCtrd.ecoreg %>% 
  datatable(filter = 'top')

AA area summary by STATE

aaCtrd.ecoreg %>%
  group_by(STATE) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  arrange(desc(sum.sqmiles)) %>% 
  datatable()

aaCtrd.ecoreg %>%
  group_by(STATE) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  filter(sum.sqmiles > 35) %>% 
  ggplot(aes(x= reorder(STATE,-sum.sqmiles), y = sum.sqmiles)) +
  geom_bar(stat='identity') +
  labs(x="",y="Sum AA area (square miles)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

AA area summary by USACE DIVISION

aaCtrd.ecoreg %>%
  group_by(DIVISION) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  arrange(desc(sum.sqmiles)) %>% 
  datatable()

USACE Division: Sum acres

# Graphs
# plot 1: sum of acres in AA by Division
pl01 <- aaCtrd.ecoreg %>% 
  group_by(DIVISION) %>% 
  summarise(sum.Acres = sum(AAidAcre)) %>% 
  ggplot(aes(reorder(x=DIVISION, sum.Acres),y = sum.Acres)) +
  geom_bar(stat='identity') +
  labs(x = "", y = "AA area (acres)", title = "Sum AA area") +
  coord_flip() +
  scale_y_log10(labels = scales::comma)
  # scale_y_log10()
# pl01
######## horizontal boxplot via ggstance
## # In ggstance, you supply aesthetics in their natural order:
# Horizontal with ggstance
# ggplot(mpg, aes(hwy, class, fill = factor(cyl))) +
#   geom_boxploth() 
# boxplot
pl02 <- aaCtrd.ecoreg %>%
  # filter()
  ggplot(aes(x=AAidAcre,y = DIVISION)) +
  geom_boxploth() +
  labs(x = "AA area (acres)", y = "", title = "Individual AA area") +
  # coord_flip() +
  scale_x_log10(labels = scales::comma)
  # scale_x_continuous(labels=fancy_scientific) ## this uses the function at the head of the doc
  # scale_x_continuous(labels=function(n){format(n, scientific = FALSE)}) # this also works but log
# pl02
## save as a grid
# plot_grid(pl01, pl02, labels = c("A", "B"))
p3 <- cowplot::plot_grid(pl01, pl02, labels = c("A", "B"), nrow = 2, align = "V")
# save_plot("v15_AA_boplots_totalarea_indi.png", p3, ncol = 2)
p3

## plot the AA area
# med.div <- aaCtrd.ecoreg %>% 
#   group_by(DIVISION, District) %>% 
#   summarise(med.AAha = median(AAareaHa), mean.areaHa = mean(AAareaHa), med.AAacre = median(AAidAcre), mean.AAacre = mean(AAidAcre))
# 
# med.div %>% 
#   ggplot(aes(x = reorder(District, med.AAacre), y = med.AAacre)) +
#   geom_point() +
#   facet_wrap(~DIVISION, scales = 'free') +
#   scale_y_log10()

AA area summary by USACE District

aaCtrd.ecoreg %>%
  group_by(District) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  arrange(desc(sum.sqmiles)) %>% 
  datatable()

USACE Division: Boxplot individual AA acres

# aaCtrd.ecoreg %>% 
#   group_by(District) %>% 
#   summarise(sum.Acres = sum(AAidAcre)) %>% 
#   ggplot(aes(reorder(x=District, sum.Acres),y = sum.Acres)) +
#   geom_bar(stat='identity') +
#   labs(x = "", y = "AA area (acres)", title = "Sum AA area") +
#   # coord_flip() +
#   scale_y_log10()
pl04 <- aaCtrd.ecoreg %>% 
  ggplot(aes(x=AAidAcre,y = District, fill= District)) +
  geom_boxploth() +
  labs(x = "Area (acres)", y = "", title = "Individual AA area by USACE District") +
  facet_wrap(~DIVISION, scales='free_y', ncol=2)+
  guides(fill=FALSE) +
  scale_x_log10(labels = scales::comma)
# if want to hide ALL legends: theme(legend.postion = "none")
pl04

# ggsave(filename = "v15_AAdistAreaBox.png",plot = pl04, width = 8, height = 8, dpi = 300)

AA area summary by PROVINCE

aaCtrd.ecoreg %>%
  group_by(PROVINCE) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  arrange(desc(sum.sqmiles)) %>% 
  datatable()

Analysis of AA NWI metadata

Image year for all AA

######
nwi.AAv15.summary <- aaCtrd.ecoreg %>%
  group_by(IMAGE_YR) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  filter(IMAGE_YR != "0") %>% 
  arrange(desc(IMAGE_YR))
datatable(nwi.AAv15.summary)

img.pl1 <- nwi.AAv15.summary %>% 
  ggplot(aes(x = IMAGE_YR, y = sum.aaAcre)) +
  geom_bar(stat='identity', color = 'black', fill = "ivory2") +
  labs(x = "NWI image year", y = "Area (acres)", title = "")
# img.pl1
img.pl2 <- nwi.AAv15.summary %>% 
  ggplot(aes(x = IMAGE_YR, y = cnt)) +
  geom_bar(stat='identity', color = 'black', fill = "ivory2") +
  labs(x = "NWI image year", y = "AA count", title = "")
# img.pl2
## save as a grid
# plot_grid(pl01, pl02, labels = c("A", "B"))
img.pl3 <- cowplot::plot_grid(img.pl1, img.pl2, labels = c("A", "B"), nrow = 1, align = "H")
img.pl3

# save_plot("v15_NWImeta_2panel.png", img.pl3, ncol = 2)

Image year by Division AA

# **Image year for all AA**
nwi.ImgYr.DIV.summary <- aaCtrd.ecoreg %>%
  group_by(IMAGE_YR, DIVISION) %>%
  summarize(cnt = n(), sum.aaAcre = sum(AAidAcre), sum.aaHa = sum(AAareaHa)) %>%
  mutate(sum.sqmiles = 0.00386102*sum.aaHa) %>% 
  filter(IMAGE_YR != "0") %>% 
  arrange(desc(IMAGE_YR))
datatable(nwi.ImgYr.DIV.summary)

imgDiv.pl1 <- nwi.ImgYr.DIV.summary %>% 
  ggplot(aes(x = IMAGE_YR, y = sum.aaAcre)) +
  geom_bar(stat='identity', color = 'black', fill = "ivory2") +
  labs(x = "NWI image year", y = "Area (acres)", title = "") +
  facet_wrap(~DIVISION, scales = 'free_y', ncol=2) +
  theme_smFacet
imgDiv.pl1

# ggsave("NWIimgYr_Divis_Acres.png", width = 6, height = 5, dpi = 300)
imgDiv.pl2 <- nwi.ImgYr.DIV.summary %>% 
  ggplot(aes(x = IMAGE_YR, y = cnt)) +
  geom_bar(stat='identity', color = 'black', fill = "ivory2") +
  labs(x = "NWI image year", y = "AA count", title = "") +
  facet_wrap(~DIVISION, scales = 'free_y', ncol=2) +
  theme_smFacet
imgDiv.pl2

# ggsave("NWIimgYr_Divis_Cnt.png", width = 6.5, height = 6, dpi = 300)
# ## save as a grid
# # plot_grid(pl01, pl02, labels = c("A", "B"))
# imgDiv.pl3 <- cowplot::plot_grid(imgDiv.pl1, imgDiv.pl2, labels = c("A", "B"), nrow = 1, align = "H")
# imgDiv.pl3
# # save_plot("v15_NWImeta_2panel.png", img.pl3, ncol = 2)

Analysis of wetland area: 3 factor and full NWI

Summary of wetland area across AA

v15NWIallType_3F_UNI <- read_csv("D:/Dropbox/PROJECTS/CORPS_NWI3P/RStudio_Projects/NWI_Analysis2016/data/Final_shapefile/V15_data/v15NWIallType_3F_UNI.csv")

Parsed with column specification:
cols(
  OBJECTID = col_integer(),
  AAv14id = col_character(),
  PolyAreaHa = col_double(),
  Shape_Length = col_double(),
  Shape_Area = col_double(),
  PolyTyUni = col_character()
)

# View(v15NWIallType_3F_UNI)
# join the aa
wetl.uni01 <- left_join(v15NWIallType_3F_UNI, aaCtrd.ecoreg, by = "AAv14id")
## quick check
# wetl.uni01 %>% 
#   names()
#  [1] "OBJECTID.x"   "AAv14id"      "PolyAreaHa"   "Shape_Length" "Shape_Area.x"
#  [6] "PolyTyUni"    "OBJECTID.y"   "Shape_Leng"   "Shape_Area.y" "AAidAcre"    
# [11] "AAareaHa"     "CenLat"       "CenLon"       "District"     "DIVISION"    
# [16] "US_L4CODE"    "US_L4NAME"    "US_L3NAME"    "NA_L3NAME"    "NA_L2NAME"   
# [21] "BlyDOMAIN"    "BlyDIVISION"  "PROVINCE"     "SECTION"      "STATE"       
# [26] "isoBioCli"
## total wetland area
# wetl.uni01 %>% 
#   names()
# unique(wetl.uni01$PolyTyUni)

Summarize area across all sites

###
uni.all <- wetl.uni01 %>% 
  # names() %>%
  group_by(PolyTyUni) %>%
  summarize(sum.areaHa = sum(PolyAreaHa)) %>%  
  mutate(Acres = sum.areaHa*2.47105) %>% 
  mutate(SqMile = Acres*0.0015625)
uni.all %>% datatable()

uni.all %>% 
  ggplot(aes(x = Acres, y = PolyTyUni)) +
  geom_barh(stat='identity') +
  labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "All mapped Cowardin types") +
  theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma)

# ggsave("AllCowTypeNWI_3F_union_Bar.png", width = 5.5, height = 4, dpi = 300)

###
AAv15.uni_wide1 <- wetl.uni01 %>% 
  group_by(AAv14id, PolyTyUni) %>%
  summarize(sum.areaHa = sum(PolyAreaHa)) %>%
  spread(key = PolyTyUni, value = sum.areaHa) %>% 
  rename(Wet3FonlyHa = Wet3Fonly) %>% 
  rename(WetNWIonlyHa = WetNWIonly) %>% 
  rename(WetUnionHa = WetUnion)
# this is only 1538 aaid. are the others with no area? need to check
# join with the AA centroid and ecoregion
AAv15.uni_wide1 <- left_join(aaCtrd.ecoreg, AAv15.uni_wide1, by = "AAv14id")
### replace NA with 0
AAv15.uni_wide1[is.na(AAv15.uni_wide1)] <- 0 # this replaces NA with 0 for all Columns
# ## SO example of changing for specific columns...
# x[c("a", "b")][is.na(x[c("a", "b")])] <- 0
## GOF mapcurves
AAv15.uni_wide1 <- AAv15.uni_wide1 %>% 
  mutate(gof = ((WetUnionHa)/(WetUnionHa + WetNWIonlyHa))*((WetUnionHa)/(WetUnionHa + Wet3FonlyHa)))
## KEY calculations: 
AAv15.uni_wide1 <- AAv15.uni_wide1 %>% 
  # names() %>% 
  mutate(TotWetHa.nwi = WetNWIonlyHa + WetUnionHa) %>%
  mutate(TotWetHa.3f = Wet3FonlyHa + WetUnionHa) %>% 
  mutate(PercWetNWI = 100*TotWetHa.nwi/AAareaHa) %>% 
  mutate(PercWet3f = 100*TotWetHa.3f/AAareaHa) %>% 
  mutate(NDAI = (TotWetHa.3f-TotWetHa.nwi)/(TotWetHa.3f+TotWetHa.nwi))
AAv15.uni_wide1 %>% 
  # names() %>% 
  datatable() %>% 
  formatRound(columns=c("TotWetHa.nwi","TotWetHa.3f","PercWetNWI","PercWet3f", "NDAI","gof"), digits=2)

Division tidy

### Calcualte acres an tidy and plot
div.Acres <- AAv15.uni_wide1 %>% 
  dplyr::select(WetUnionHa, Wet3FonlyHa, WetNWIonlyHa, DIVISION) %>% 
  rename(Union = WetUnionHa, Only3F = Wet3FonlyHa, OnlyNWI = WetNWIonlyHa) %>% 
  gather(key = uniType, value = AreaHa, -DIVISION) %>% 
  mutate(AreaAcre = AreaHa*2.47105) %>% 
  group_by(uniType, DIVISION) %>%
  summarise(sum.Acre = sum(AreaAcre))
datatable(div.Acres)

###
div.Acres %>% 
  ggplot(aes(x= sum.Acre, y = uniType)) +
  geom_barh(stat='identity') +
  facet_wrap(~DIVISION, scales = "free") +
  labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "All mapped Cowardin types") +
  theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma)

div.Acres %>% 
  ggplot(aes(x= sum.Acre, y = DIVISION, fill = uniType)) +
  geom_barh(stat='identity', color = 'black') +
  facet_wrap(~DIVISION, scales = "free", ncol = 1) +
  labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "All mapped Cowardin types") +
  theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma) +
  theme(axis.ticks = element_blank(), axis.text.y = element_blank())

# ggsave("DivisionUniType.png", width = 6.5, height = 6, dpi = 300)
# uni.all %>% 
#   ggplot(aes(x = Acres, y = PolyTyUni)) +
#   geom_barh(stat='identity') +
#   labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "All mapped Cowardin types") +
#   theme_bw() +
#   # scale_x_log10(labels = scales::comma)
#   scale_x_continuous(labels = scales::comma)

### Calcualte acres an tidy and plot
## Baily provinces
Prov.Acres <- AAv15.uni_wide1 %>% 
  dplyr::select(WetUnionHa, Wet3FonlyHa, WetNWIonlyHa, PROVINCE) %>% 
  rename(Union = WetUnionHa, Only3F = Wet3FonlyHa, OnlyNWI = WetNWIonlyHa) %>% 
  gather(key = uniType, value = AreaHa, -PROVINCE) %>% 
  mutate(AreaAcre = AreaHa*2.47105) %>% 
  group_by(uniType, PROVINCE) %>%
  summarise(sum.Acre = sum(AreaAcre))
datatable(Prov.Acres)

###
Prov.Acres %>% 
  ggplot(aes(x= sum.Acre, y = uniType)) +
  geom_barh(stat='identity') +
  facet_wrap(~PROVINCE, scales = "free") +
  labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "All mapped Cowardin types") +
  theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma)

### filter out sites >50 acreas in AA size
Prov.Acres %>% 
  group_by(PROVINCE) %>% 
  mutate(sumAllPRovinceArea = sum(sum.Acre)) %>% 
  filter(sumAllPRovinceArea > 50) %>% 
  ggplot(aes(x= sum.Acre, y = PROVINCE, fill = uniType)) +
  geom_barh(stat='identity', color = 'black') +
  facet_wrap(~PROVINCE, scales = "free", ncol = 2) +
  labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "All mapped Cowardin types", fill="") +
  theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma) +
  theme(axis.ticks = element_blank(), axis.text.y = element_blank()) 

# ggsave("PROVINCWUniType.png", width = 12.5, height = 10, dpi = 300)

Percent of AA as NWI (all C types) and 3 factor

Summary: percent AA as NWI

## NWI
summary(AAv15.uni_wide1$PercWetNWI)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   2.911  14.020  17.280 100.000 

# dir()
# AAv15.uni_wide1 %>% 
#   summarise_at(.cols = PercWet3f,.funs = summary)

Summary: percent AA as 3F

summary(AAv15.uni_wide1$PercWet3f)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.0000   0.1669   3.0570  12.8500  13.6100 100.0000 

## all AA
AAv15.uni_wide1 %>% 
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.05, color = "blue") +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "All NWI Cowardin types") +
  theme_bw()

  # geom_hex()
# ggsave("AllNWICowardinTy_AllAA_percentAAWet.png", width = 5, height = 5, dpi = 300)

% AA full NWI by DIVISION

## AA by Division full NWI
AAv15.uni_wide1 %>% 
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.05, color = "blue") +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "Full NWI") +
  facet_wrap(~DIVISION, ncol = 4) +
  xlim(0,100) +
  ylim(0,100) +
  theme(strip.text.x = element_text(size = 8))

      
# ggsave(filename = "PercentAA_NWIallType.png", width = 5, height = 10, dpi = 300)
# ggsave(filename = "PercentAA_NWIallTypeW.png", width = 10, height = 5, dpi = 300)

AAv15.uni_wide1 %>%
  group_by(IMAGE_YR) %>% 
  summarize(mean.PercWetNWI = mean(PercWetNWI), sd.PercWetNWI = sd(PercWetNWI), mean.PercWet3f = mean(PercWet3f), sd.PercWet3f = sd(PercWet3f)) %>% 
  filter(IMAGE_YR != 0) %>% 
  mutate(YrBin = ifelse(IMAGE_YR <= 1979, "1970's",
               ifelse(IMAGE_YR >= 1980 & IMAGE_YR < 1990, "1980's",ifelse(IMAGE_YR >= 1990 & IMAGE_YR <2000, "1990's", "2000's")))) %>%
  ggplot(aes(x = mean.PercWetNWI, y = mean.PercWet3f)) + 
  geom_point(alpha = 0.25, color = "blue", size = 6) +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "") +
  facet_wrap(~YrBin, ncol = 2, scales = "free") +
  # xlim(0,100) +
  # ylim(0,100) +
  theme(strip.text.x = element_text(size = 8)) +
  theme_smFacet

####
AAv15.uni_wide1 %>%
  # group_by(IMAGE_YR) %>% 
  # summarize(mean.PercWetNWI = mean(PercWetNWI), sd.PercWetNWI = sd(PercWetNWI), mean.PercWet3f = mean(PercWet3f), sd.PercWet3f = sd(PercWet3f)) %>% 
  filter(IMAGE_YR != 0) %>% 
  mutate(YrBin = ifelse(IMAGE_YR <= 1979, "1970's",
               ifelse(IMAGE_YR >= 1980 & IMAGE_YR < 1990, "1980's",ifelse(IMAGE_YR >= 1990 & IMAGE_YR <2000, "1990's", "2000's")))) %>%
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.085, color = "blue", size = 3) +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "") +
  facet_wrap(~YrBin, ncol = 2, scales = "free") +
  # xlim(0,100) +
  # ylim(0,100) +
  theme(strip.text.x = element_text(size = 8)) +
  theme_smFacet

# ggsave("NWI_ImgBin_scatterPercent.png", width = 6, height = 5, dpi = 300)

Percent of AA as NWI (REDUCED C types) and 3 factor

This ihas been rduced by eliminating NWI polygons in the followong Cowardin classes: 1. Riverine 1. Lacustrine 1. Other 1. Marine Deep Water

t1 <- gsub("\\.shp$","", list.files("data/Final_shapefile/v15_data", pattern="shp$", full.names=TRUE, recursive = FALSE))
# t2 <- gsub("\\.shp$","", dir("data/Final_shapefile/v15_data", pattern="shp$", full.names=TRUE, recursive = FALSE)) # Note list.files and dir seem to do the exact smae thing...
nwiRed01.tbl <- map2(.x = "data/Final_shapefile/v15_data", .y = "FinNWI_noPoLacRiv_3Funi", .f =  readOGR) 

OGR data source with driver: ESRI Shapefile 
Source: "data/Final_shapefile/v15_data", layer: "FinNWI_noPoLacRiv_3Funi"
with 49274 features
It has 14 fields

Z-dimension discarded

nwiRed01.tbl <- as_tibble(nwiRed01.tbl[[1]]@data) 
names(nwiRed01.tbl)

 [1] "OBJECTID"   "FID_FinNWI" "AAv14id"    "WETLAND_TY" "PolyAreaM2"
 [6] "PolyAreaHa" "PolyType"   "FID_Fin3F_" "AAv14id_1"  "PolyType_1"
[11] "Shape_Leng" "Shape_Area" "PolyTyUni"  "UniClass"  

# aacentr.tbl %>% datatable()
## summarize total AA area
nwiRed01.tbl %>% 
  # names() %>%
  summarize(sum.areaHa = sum(PolyAreaHa), cnt.AApoly = n()) %>%
  mutate(sum.Acres = 2.47105*sum.areaHa) %>% 
  mutate(sum.sqmiles = 0.00386102*sum.areaHa) %>% 
  datatable()

# join the aa, eroregion stuff
wetl.uni01red <- left_join(nwiRed01.tbl, aaCtrd.ecoreg, by = "AAv14id")

joining character vector and factor, coercing into character vector

## reduced NWI
uni.all.red <- wetl.uni01red %>% 
  # names() %>%
  group_by(PolyTyUni) %>%
  summarize(sum.areaHa = sum(PolyAreaHa)) %>%  
  mutate(Acres = sum.areaHa*2.47105) %>% 
  mutate(SqMile = Acres*0.0015625)
uni.all.red %>% datatable()

uni.all.red  %>% 
  ggplot(aes(x = Acres, y = PolyTyUni)) +
  geom_barh(stat='identity') +
  labs(y = "", x = "Area (acres)", fill = "", title = "Mapped area by geometric union", subtitle = "Reduced set of Cowardin types") +
  # theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma) +
  theme_CowsMod1

wetl.uni01red %>% 
  # names() %>%
  group_by(PolyTyUni) %>%
  summarize(sum.areaHa = sum(PolyAreaHa)) %>%  
  mutate(WetAcres = sum.areaHa*2.47105) %>% 
  mutate(WetSqMile = WetAcres*0.0015625)

DIVISION: NWI reduced Acreas stacked bar

### reduced wetland type
## DIVSION sum 
red.tdy <- left_join(aaCtrd.ecoreg, nwiRed01.tbl, by = "AAv14id")

joining factor and character vector, coercing into character vector

red.tdy %>% 
  group_by(PolyTyUni) %>% 
  summarise(sum.areaHa = sum(PolyAreaHa)) %>% 
  mutate(sum.acre = sum.areaHa*2.47105)

div.Acres.red <- red.tdy %>% 
  dplyr::select(PolyTyUni, DIVISION, PolyAreaHa) %>% 
  mutate(AreaAcre = PolyAreaHa*2.47105) %>% 
  na.omit() %>% 
  group_by(PolyTyUni, DIVISION) %>%
  summarise(sum.Acre = sum(AreaAcre))
# datatable(div.Acres)
###
div.Acres.red %>%
  ggplot(aes(x= sum.Acre, y = PolyTyUni)) +
  geom_barh(stat='identity') +
  facet_wrap(~DIVISION, scales = "free", ncol = 4)

div.Acres.red %>% 
  ggplot(aes(x= sum.Acre, y = DIVISION, fill = PolyTyUni)) +
  geom_barh(stat='identity', color = 'black') +
  facet_wrap(~DIVISION, scales = "free", ncol = 1) +
  labs(y = "", x = "Area (acres)",fill = "") +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma) +
  theme(axis.ticks = element_blank(), axis.text.y = element_blank())

# ggsave("DivisionUniType_Reduced.png", width = 6.5, height = 6, dpi = 300)

# reduced nwi
wetl.uni01red_wide1 <- wetl.uni01red %>% 
  group_by(AAv14id, PolyTyUni) %>%
  summarize(sum.areaHa = sum(PolyAreaHa)) %>%
  spread(key = PolyTyUni, value = sum.areaHa) %>% 
  rename(Wet3FonlyHa = Wet3Fonly) %>% 
  rename(WetNWIonlyHa = WetNWIonly) %>% 
  rename(WetUnionHa = WetUnion)
# this is only 1499 (from 1538 aaid for the nwi full 3f uni layer). 
# Need to check
# join with the AA centroid and ecoregion
wetl.uni01red_wide1 <- left_join(aaCtrd.ecoreg, wetl.uni01red_wide1, by = "AAv14id")
### replace NA with 0
wetl.uni01red_wide1[is.na(wetl.uni01red_wide1)] <- 0 # this replaces NA with 0 for all Columns
# ## SO example of changing for specific columns...
# x[c("a", "b")][is.na(x[c("a", "b")])] <- 0
## GOF mapcurves
wetl.uni01red_wide1 <- wetl.uni01red_wide1 %>% 
  mutate(gof = ((WetUnionHa)/(WetUnionHa + WetNWIonlyHa))*((WetUnionHa)/(WetUnionHa + Wet3FonlyHa)))
## KEY calculations: 
wetl.uni01red_wide1 <- wetl.uni01red_wide1 %>% 
  # names() %>% 
  mutate(TotWetHa.nwi = WetNWIonlyHa + WetUnionHa) %>%
  mutate(TotWetHa.3f = Wet3FonlyHa + WetUnionHa) %>% 
  mutate(PercWetNWI = 100*TotWetHa.nwi/AAareaHa) %>% 
  mutate(PercWet3f = 100*TotWetHa.3f/AAareaHa) %>% 
  mutate(NDAI = (TotWetHa.3f-TotWetHa.nwi)/(TotWetHa.3f+TotWetHa.nwi))
wetl.uni01red_wide1 %>% 
  # names() %>% 
  datatable() %>% 
  formatRound(columns=c("TotWetHa.nwi","TotWetHa.3f","PercWetNWI","PercWet3f", "NDAI","gof"), digits=2)

## all AA: comparison of reduced and full nwi
b1 <- AAv15.uni_wide1 %>% 
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.05, color = "blue") +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "All NWI Cowardin types") +
  theme_bw()
  # geom_hex()
b2 <- wetl.uni01red_wide1 %>% 
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.05, color = "blue") +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "Reduced NWI Cowardin types") +
  theme_bw()
  # geom_hex()
# ggsave("ReducedNWICowardinTy_AllAA_percentAAWet.png", width = 5, height = 5, dpi = 300)
b3 <- plot_grid(b1, b2, labels = c("A", "B"), align = "h")
# save_plot("v15_NWIred.png", b3, ncol = 2)

Reduced: difference in %AA

• Summary of the %AA for all sites*

## calc the diff
wetl.uni01red_wide1 <- wetl.uni01red_wide1 %>% 
  mutate(diff.AApercWet = PercWet3f-PercWetNWI)
Some.wet.aa <- wetl.uni01red_wide1 %>%
  mutate(combined_WetArea = TotWetHa.3f + TotWetHa.nwi) %>% 
  filter(combined_WetArea >0)
summary(Some.wet.aa$diff.AApercWet)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-100.0000   -3.4770    0.7392    1.6750    6.5000  100.0000 

##
areaZ.wide <- wetl.uni01red_wide1 %>% 
  dplyr::select(AAv14id, AAareaHa, PercWet3f, PercWetNWI,TotWetHa.nwi, TotWetHa.3f)
areaZ.tidy <- areaZ.wide %>% 
  gather(key = PropVar, value = Prop, PercWet3f:PercWetNWI) %>% 
  gather(key = AreaHaVar, value = AreaHa, TotWetHa.nwi:TotWetHa.3f)
areaZ.tidy %>% 
  ggplot(aes(x=AreaHaVar, y = AreaHa)) +
  geom_boxplot()

areaZ.tidy %>% 
  ggplot(aes(x=PropVar, y = Prop)) +
  geom_boxplot()

## AA by Division: REDUCED
wetl.uni01red_wide1 %>% 
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.05, color = "blue") +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "Proportion of AA mapped by 3F and NWI (reduced Cowardin types)") +
  facet_wrap(~DIVISION, ncol = 4) +
  xlim(0,100) +
  ylim(0,100) +
  # theme(strip.text.x = element_text(size = 8))
  theme_smFacet

  # theme_CowsMod1
      
# ggsave(filename = "PercentAA_NWIreducedTypeN.png", width = 5, height = 10, dpi = 300)
# ggsave(filename = "PercentAA_NWIreducedTypeW.png", width = 10, height = 5, dpi = 300)

wetl.uni01red_wide1 %>%
  group_by(IMAGE_YR) %>% 
  summarize(mean.PercWetNWI = mean(PercWetNWI), sd.PercWetNWI = sd(PercWetNWI), mean.PercWet3f = mean(PercWet3f), sd.PercWet3f = sd(PercWet3f)) %>% 
  filter(IMAGE_YR != 0) %>% 
  mutate(YrBin = ifelse(IMAGE_YR <= 1979, "1970's",
               ifelse(IMAGE_YR >= 1980 & IMAGE_YR < 1990, "1980's",ifelse(IMAGE_YR >= 1990 & IMAGE_YR <2000, "1990's", "2000's")))) %>%
  ggplot(aes(x = mean.PercWetNWI, y = mean.PercWet3f)) + 
  geom_point(alpha = 0.25, color = "blue", size = 6) +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "") +
  facet_wrap(~YrBin, ncol = 2, scales = "free") +
  # xlim(0,100) +
  # ylim(0,100) +
  theme(strip.text.x = element_text(size = 8)) +
  theme_smFacet

####
wetl.uni01red_wide1 %>%
  # group_by(IMAGE_YR) %>% 
  # summarize(mean.PercWetNWI = mean(PercWetNWI), sd.PercWetNWI = sd(PercWetNWI), mean.PercWet3f = mean(PercWet3f), sd.PercWet3f = sd(PercWet3f)) %>% 
  filter(IMAGE_YR != 0) %>% 
  mutate(YrBin = ifelse(IMAGE_YR <= 1979, "1970's",
               ifelse(IMAGE_YR >= 1980 & IMAGE_YR < 1990, "1980's",ifelse(IMAGE_YR >= 1990 & IMAGE_YR <2000, "1990's", "2000's")))) %>%
  ggplot(aes(x = PercWetNWI, y = PercWet3f)) + 
  geom_point(alpha = 0.085, color = "blue", size = 3) +
  geom_smooth(color = "red", method = "lm") +
  geom_abline(slope = 1, intercept = 0, lty = 'dashed', size = 1.1, alpha = .4) +
  # geom_density_2d() +
  labs(x = "% AA mapped as NWI", y = "% AA mapped as 3 factor", title = "Reduced Cowardin types") +
  facet_wrap(~YrBin, ncol = 2, scales = "free") +
  # xlim(0,100) +
  # ylim(0,100) +
  theme(strip.text.x = element_text(size = 8)) +
  theme_smFacet

# ggsave("NWI_ImgBin_scatterPercent_reduced.png", width = 6, height = 5, dpi = 300)

Division: NDAI reduced NWI

pl06 <- wetl.uni01red_wide1 %>%
  group_by(DIVISION) %>% 
  mutate(mean.NDAI = mean(NDAI, na.rm = TRUE)) %>% 
  # names() %>% 
  ggplot(aes(x=NDAI, y = reorder(DIVISION, mean.NDAI))) +
  geom_boxploth(fill = "ivory2") +
  # geom_jitter(color = "blue", alpha=0.05, size = 2) +
  geom_vline(xintercept = 0, size=1.5, color = "red", lty = "dashed", alpha = 0.5) +
  labs(y="", title = "Normalized Difference Area Index", subtitle = "Reduced set of Cowardin types") +
  theme_CowsMod1
# # save_plot(plot = pl06, filename = "DIVISION_NDAI_boxplot_Red.png")
# ggsave(filename = "Division_NDAI_boxplotNoJit_red.png", plot = pl06, dpi = 300, width = 7.5, height = 4)
pl06

District: NDAI

pl07 <- wetl.uni01red_wide1 %>% 
  # mutate(mean.NDAI = mean(NDAI, na.rm = TRUE)) %>%
  group_by(District) %>% 
  mutate(med.NDAI = median(NDAI, na.rm = TRUE)) %>%
  mutate(cnt = n()) %>% 
  filter(cnt >= 15) %>% 
  # names() %>% 
  ggplot(aes(x=NDAI, y = reorder(District,med.NDAI))) +
  geom_boxploth(fill = "ivory2") +
  # geom_jitter(color = "blue", alpha=0.05, size = 2) +
  geom_vline(xintercept = 0, size=1.5, color = "red", lty = "dashed", alpha = 0.5) +
  labs(y="", title = "Normalized Difference Area Index", subtitle = "Reduced set of Cowardin types") +
  theme_CowsMod1
# ggsave(filename = "District_NDAI_boxplotNojitter_red.png", plot = pl07, dpi = 300, width = 7.5, height = 6)
pl07

Province: NDAI (#fail)

## can't get this shit to work...

pl08 <- wetl.uni01red_wide1 %>% 
  # mutate(mean.NDAI = mean(NDAI, na.rm = TRUE)) %>%
  group_by(PROVINCE) %>% 
  mutate(med.NDAI = median(NDAI, na.rm = TRUE)) %>%
  mutate(cnt = n()) %>% 
  filter(cnt >= 15) %>% 
  # names() %>% 
  ggplot(aes(x=NDAI, y = reorder(PROVINCE,med.NDAI))) +
  geom_boxploth(fill = "ivory2") +
  # geom_jitter(color = "blue", alpha=0.05, size = 2) +
  geom_vline(xintercept = 0, size=1.5, color = "red", lty = "dashed", alpha = 0.5) +
  labs(y="", title = "Normalized Difference Area Index", subtitle = "Reduced set of Cowardin types") +
  theme_smFacet +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) +
  # theme_CowsMod1
  facet_wrap(~PROVINCE, ncol = 1)


ggsave(filename = "BlyProvince_NDAI_boxplot_red.png", plot = pl08, dpi = 300, width = 8.5, height = 15)

pl08

####
disti.ndai.sum <- wetl.uni01red_wide1 %>% 
  group_by(District) %>% 
  summarise(cnt = n(), mean.ndai = mean(NDAI, na.rm=TRUE), med.ndai = median(NDAI, na.rm=TRUE), q25.ndai = quantile(NDAI,prob = .25, na.rm=TRUE), q75.ndai = quantile(NDAI,prob = .75, na.rm=TRUE))
disti.ndai.sum %>% 
  datatable() %>% 
  formatRound(~mean.ndai + med.ndai + q25.ndai + q75.ndai, 2)

## error bar with median NDAI
disti.ndai.sum %>% 
  filter(cnt >= 30) %>%
  ggplot(aes(x=med.ndai, y = reorder(District, med.ndai))) +
  geom_errorbarh(aes(xmin= q25.ndai, xmax = q75.ndai)) +
  geom_point(size = 4.5, color = "red")

########
## Note cool use of summarize_at
########
wetl.uni01red_wide1 %>% 
  summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd', 'median'))

## DIVISION
wetl.uni01red_wide1 %>% 
  group_by(DIVISION) %>%
  summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd')) %>% 
  datatable() %>% 
  formatRound(~AAareaHa_mean + PercWet3f_mean + PercWetNWI_mean + AAareaHa_sd + PercWet3f_sd + PercWetNWI_sd, 2)

## DISTRICT
wetl.uni01red_wide1 %>% 
  group_by(District) %>%
  summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd')) %>% 
  # names() %>% 
  datatable() %>% 
  formatRound(~AAareaHa_mean + PercWet3f_mean + PercWetNWI_mean + AAareaHa_sd + PercWet3f_sd + PercWetNWI_sd, 2)

  
## US_L3NAME
wetl.uni01red_wide1 %>% 
  group_by(US_L3NAME) %>%
  summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd')) %>% 
  # names() %>% 
  datatable() %>% 
  formatRound(~AAareaHa_mean + PercWet3f_mean + PercWetNWI_mean + AAareaHa_sd + PercWet3f_sd + PercWetNWI_sd, 3)

Difference in the %AA mapped as wetland

## historgram of the $ diff in % AA mapped
# wetl.uni01red_wide1 %>% 
#   ggplot(aes(x=diff.AApercWet)) +
#   geom_histogram()
wetl.uni01red_wide1 %>% 
  ggplot(aes(x=diff.AApercWet)) +
  geom_density(fill = 'ivory2') +
  geom_vline(xintercept = 0, size = 1.2, lty = "dashed", color = "red", alpha = .2) +
  labs(x = "Difference in %AA mapped", y = "Density")

# ggsave("DiffInPercAAasWet.png", width = 5.5, height = 5, dpi = 300)
## DIVISION
wetl.uni01red_wide1 %>% 
  group_by(DIVISION) %>%
  summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd')) %>% 
  ggplot(aes(x = diff.AApercWet_mean
, y = reorder(DIVISION,diff.AApercWet_mean))) +
geom_point()

## plot faceted version of graph
wetl.uni01red_wide1 %>% 
  ggplot(aes(x=diff.AApercWet)) +
  geom_density(fill = 'ivory2') +
  geom_vline(xintercept = 0, size = 1.2, lty = "dashed", color = "red", alpha = .2) +
  labs(x = "Difference in %AA mapped", y = "Density") +
  facet_wrap(~DIVISION, ncol = 2, scales = 'free_y') +
  theme_smFacet

# ggsave("DiffInPercAAasWet_DIVden.png", width = 7.5, height = 6, dpi = 300)
# ## DISTRICT
# wetl.uni01red_wide1 %>% 
#   group_by(District) %>%
#   summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd')) %>% 
#   # names() %>% 
#   datatable() %>% 
#   formatRound(~AAareaHa_mean + PercWet3f_mean + PercWetNWI_mean + AAareaHa_sd + PercWet3f_sd + PercWetNWI_sd, 2)
#   
# 
# ## US_L3NAME
# wetl.uni01red_wide1 %>% 
#   group_by(US_L3NAME) %>%
#   summarise_at(c('AAareaHa','PercWet3f','PercWetNWI','diff.AApercWet'),c('mean','sd')) %>% 
#   # names() %>% 
#   datatable() %>% 
#   formatRound(~AAareaHa_mean + PercWet3f_mean + PercWetNWI_mean + AAareaHa_sd + PercWet3f_sd + PercWetNWI_sd, 3)

## NDAI 
## plot: reduced NWI
wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  ggplot(aes(x=NDAI)) +
  geom_density(alpha = 0.4, fill = 'ivory2', color = "black") +
  geom_vline(xintercept = 0, color = 'red', size = 1.15, lty = "dashed") +
  facet_wrap(~DIVISION, ncol = 2, scale = 'free_y') +
  labs(x="NDAI", y = "Density", title = "NDAI by USACE Division") +
  theme(legend.position = "none")

# ggsave(file="NDAI_Division_Density_red.png", width = 7.5, height = 6, dpi = 300)
## histogram: reduced NWI
wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  ggplot(aes(x=NDAI)) +
  geom_histogram(alpha = 0.4, fill = 'ivory2', color = "black") +
  geom_vline(xintercept = 0, color = 'red', size = 1.15, lty = "dashed") +
  facet_wrap(~DIVISION, ncol = 2, scale = 'free_y') +
  labs(x="NDAI", y = "AA count", title = "NDAI by USACE Division") +
  theme(legend.position = "none")

# ggsave(file="NDAI_Division_histogram_red.png", width = 7.5, height = 6, dpi = 300)

USACE Division – Mapcurves GOF

gof.div.plot <- wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  ggplot(aes(x=gof)) +
  geom_density(color = "black", alpha = 0.45, fill = "ivory2") +
  # geom_vline(xintercept = 0, color = 'red', size = 1.15, lty = "dashed") +
  facet_wrap(~DIVISION, ncol = 2, scale = 'free') +
  labs(x="Goodness of fit", y = "Density", title = "Mapcurves goodness of fit by USACE Division") +
  theme(legend.position = "none")
# gof.div.plot %>% 
#   ggsave(file = "GOF_Division_red.png", width = 7, height = 7, dpi = 300)
## histgram
gof.div.hist <- wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  ggplot(aes(x=gof)) +
  geom_histogram(color = "black", alpha = 0.45, fill = "ivory2") +
  # geom_vline(xintercept = 0, color = 'red', size = 1.15, lty = "dashed") +
  facet_wrap(~DIVISION, ncol = 2, scale = 'free_y') +
  labs(x="Goodness of fit", y = "AA count", title = "Mapcurves goodness of fit by USACE Division") +
  theme(legend.position = "none")
# ggsave(file="GOF_Division_histogram_red2.png", width = 6.85, height = 6, dpi = 300)
gof.div.hist

## boxplots
wetl.uni01red_wide1 %>% 
  geom_histogram(color = "black", alpha = 0.45, fill = "ivory2") +
  # geom_vline(xintercept = 0, color = 'red', size = 1.15, lty = "dashed") +
  facet_wrap(~DIVISION, ncol = 2, scale = 'free_y') +
  labs(x="Goodness of fit", y = "AA count", title = "Mapcurves goodness of fit by USACE Division") +
  theme(legend.position = "none") +
  coord_flip()

# ggsave(file="GOF_Division_histogram_redalt.png", width = 7.5, height = 6, dpi = 300)
gof.div.hist

Baily province

Prov.Acres.red <- wetl.uni01red_wide1  %>% 
  dplyr::select(WetUnionHa, Wet3FonlyHa, WetNWIonlyHa, PROVINCE) %>% 
  rename(Union = WetUnionHa, Only3F = Wet3FonlyHa, OnlyNWI = WetNWIonlyHa) %>% 
  gather(key = uniType, value = AreaHa, -PROVINCE) %>% 
  mutate(AreaAcre = AreaHa*2.47105) %>% 
  group_by(uniType, PROVINCE) %>%
  summarise(sum.Acre = sum(AreaAcre))
datatable(Prov.Acres.red)

###
# Prov.Acres.red %>% 
#   ggplot(aes(x= sum.Acre, y = uniType)) +
#   geom_barh(stat='identity') +
#   facet_wrap(~PROVINCE, scales = "free") +
#   labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "Subset Cowardin types") +
#   theme_bw() +
#   # scale_x_log10(labels = scales::comma)
#   scale_x_continuous(labels = scales::comma)
### filter out sites >50 acreas in AA size
Prov.Acres.red %>% 
  group_by(PROVINCE) %>% 
  mutate(sumAllPRovinceArea = sum(sum.Acre)) %>% 
  filter(sumAllPRovinceArea > 50) %>% 
  ggplot(aes(x= sum.Acre, y = PROVINCE, fill = uniType)) +
  geom_barh(stat='identity', color = 'black') +
  facet_wrap(~PROVINCE, scales = "free", ncol = 2) +
  labs(y = "", x = "Area (acres)", title = "Mapped area by geometric union", subtitle = "Subset of Cowardin types", fill="") +
  theme_bw() +
  # scale_x_log10(labels = scales::comma)
  scale_x_continuous(labels = scales::comma) +
  theme(axis.ticks = element_blank(), axis.text.y = element_blank()) 

# ggsave("DivisionUniType.png", width = 9.5, height = 10, dpi = 300)

Bailey Section – Mapcurves GOF

gof.BlyProvince.plot <- wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(PROVINCE) %>% 
  mutate(cnt_by_prov  = n()) %>% 
  filter(cnt_by_prov >= 30) %>% 
  ggplot(aes(x=gof)) +
  geom_density(color = "black", alpha = 0.45, fill = "ivory2") +
  # geom_vline(xintercept = 0, color = 'red', size = 1.15, lty = "dashed") +
  facet_wrap(~PROVINCE, ncol = 1) +
  labs(x="Goodness of fit", y = "Density", title = "Mapcurves goodness of fit by Bailey Province") +
  theme(legend.position = "none") +
  theme_CowsMod1
gof.BlyProvince.plot
gof.BlyProvince.plot %>%
  ggsave(file = "GOF_BlyProvince_red.png", width = 7, height = 10, dpi = 300)

## Datatble
wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(PROVINCE) %>% 
  mutate(cnt_by_prov  = n()) %>% 
  datatable(caption = "GOF reduced NWI Baily Province")

Difference in % mapped by NWI

wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(IMAGE_YR) %>% 
  summarize(cnt_by_prov  = n(), mean.PercWetNWI = mean(PercWetNWI), mean.PercWet3f = mean(PercWet3f),sum.AAv15Acre = sum(AAidAcre)) %>% 
  mutate(dif.PercAA3F_NWI = mean.PercWet3f-mean.PercWetNWI) %>% 
  mutate(wt.mn = map2_dbl(.x = dif.PercAA3F_NWI, .y = sum.AAv15Acre, .f=weighted.mean)) %>% 
  filter(IMAGE_YR >0) %>% 
  ggplot(aes(x=IMAGE_YR, y = dif.PercAA3F_NWI)) +
  geom_hline(yintercept = 0, color='red', size=1.2, alpha=0.5) +
  geom_point(aes(size=sum.AAv15Acre), color = "blue", alpha = .3) +
  labs(x = "NWI image year", y = "Difference in % AA mapped", title = "All AA - Reduced NWI") +
  # theme_CowsMod1
  theme_smFacet +
  theme(legend.position="none") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(-100,100) 
  # geom_smooth(method = "lm", se = FALSE, lty = "dashed", color = 'black')
ggsave("DiffernceMappedNWIred_AllAA_reduced.png", width = 5.5, height = 3.5, dpi = 300)

## datatable
wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(IMAGE_YR) %>% 
  summarize(cnt_by_prov  = n(), mean.PercWetNWI = mean(PercWetNWI), mean.PercWet3f = mean(PercWet3f),sum.AAv15Acre = sum(AAidAcre)) %>% 
  mutate(dif.PercAA3F_NWI = mean.PercWet3f-mean.PercWetNWI) %>% 
  filter(IMAGE_YR >0) %>% 
  datatable(filter = 'top',caption = "mean_AAdiff")

Difference in percent AA mapped by 3f and reduced NWI DISTRICT

wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(District, IMAGE_YR) %>% 
  summarize(cnt_by_prov  = n(), mean.PercWetNWI = mean(PercWetNWI), mean.PercWet3f = mean(PercWet3f),sum.AAv15Acre = sum(AAidAcre)) %>% 
  mutate(dif.PercAA3F_NWI = mean.PercWet3f-mean.PercWetNWI) %>% 
  mutate(wt.mn = map2_dbl(.x = dif.PercAA3F_NWI, .y = sum.AAv15Acre, .f=weighted.mean)) %>% 
  filter(IMAGE_YR >0) %>% 
  ggplot(aes(x=IMAGE_YR, y = dif.PercAA3F_NWI)) +
  geom_hline(yintercept = 0, color='red', size=1.2, alpha=0.5) +
  geom_point(aes(size=sum.AAv15Acre), color = "blue", alpha = .3) +
  labs(x = "NWI image year", y = "Difference in % AA mapped", title = "Reduced NWI") +
  facet_wrap(~District) +
  # theme_CowsMod1
  theme_smFacet +
  theme(legend.position="none") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(-100,100)

# ggsave("DiffernceMappedNWIred_District.png", width = 7.5, height = 5.5, dpi = 300)

Difference in percent AA mapped by 3f and reduced NWI DIVISION

wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(DIVISION, IMAGE_YR) %>% 
  summarize(cnt_by_prov  = n(), mean.PercWetNWI = mean(PercWetNWI), mean.PercWet3f = mean(PercWet3f),sum.AAv15Acre = sum(AAidAcre)) %>% 
  mutate(dif.PercAA3F_NWI = mean.PercWet3f-mean.PercWetNWI) %>% 
  mutate(wt.mn = map2_dbl(.x = dif.PercAA3F_NWI, .y = sum.AAv15Acre, .f=weighted.mean)) %>% 
  filter(IMAGE_YR >0) %>% 
  ggplot(aes(x=IMAGE_YR, y = dif.PercAA3F_NWI)) +
  geom_hline(yintercept = 0, color='red', size=1.2, alpha=0.5) +
  geom_point(aes(size=sum.AAv15Acre), color = "blue", alpha = .3) +
  labs(x = "NWI image year", y = "Difference in % AA mapped", title = "Reduced NWI") +
  facet_wrap(~DIVISION, ncol = 4) +
  # theme_CowsMod1
  theme_smFacet +
  theme(legend.position="none") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(-100,100) +
  geom_smooth(method = "lm", se = FALSE, lty = "dashed", size = 1, color = 'black')

# ggsave("DiffernceMappedNWIred_Division.png", width = 8, height = 3.75, dpi = 300)
wetl.uni01red_wide1 %>% 
  # filter(cnt >= 7) %>% 
  group_by(DIVISION, IMAGE_YR) %>% 
  summarize(cnt_by_prov  = n(), mean.PercWetNWI = mean(PercWetNWI), mean.PercWet3f = mean(PercWet3f),sum.AAv15Acre = sum(AAidAcre)) %>% 
  mutate(dif.PercAA3F_NWI = mean.PercWet3f-mean.PercWetNWI) %>% 
  mutate(wt.mn = map2_dbl(.x = dif.PercAA3F_NWI, .y = sum.AAv15Acre, .f=weighted.mean)) %>% 
  filter(IMAGE_YR >0) %>% 
  datatable(filter = 'top')

# export a version of the wide AA with all the stuff for jpining in Arc
# wetl.uni01red_wide1 %>% 
#   write_csv("Wide_NWIred_3F_uni_CALCs_for_AAv15join.csv")

Regression modeling

Reduced NWI by DIVISION

uni.red.DIV.nDF <- wetl.uni01red_wide1 %>%
  group_by(DIVISION) %>% 
  nest()
  
uni.red.DIV.nDF <- uni.red.DIV.nDF %>% 
  mutate(model01 = purrr::map(data, ~ lm(TotWetHa.nwi ~ TotWetHa.3f, data = .))
)  
# Extract model summaries:
uni.red.DIV.nDF <- uni.red.DIV.nDF %>% unnest(model01 %>% purrr::map(broom::glance))
uni.red.DIV.nDF %>% 
  ggplot(aes(x = adj.r.squared, y= reorder(DIVISION,adj.r.squared))) +
  # geom_barh(stat = 'identity', aes(fill = p.value)) +
  geom_barh(stat = 'identity', fill = 'ivory2', color = 'ivory3') +
  labs(x=expression('adjusted R'^{2}), y = "", title = "Linear models: 3F area ~ NWI area") +
  theme_CowsMod1

# ggsave("LinerModels_Division_red.png", dpi = 300, width = 7.5, height = 3)
# # Extract model summaries:
# by_country %>% unnest(model %>% purrr::map(broom::glance))
## data table
uni.red.DIV.nDF %>% 
  datatable(filter='top')

NWI image year: adjust R2 lm

uni.red.ImgYr.nDF <- wetl.uni01red_wide1 %>%
  filter(IMAGE_YR >0) %>% 
  group_by(IMAGE_YR) %>% 
  mutate(cnt = n()) %>%
  filter(cnt>=25) %>% # note need to filter to have >1 record. I'm going for minimum n of 25 AA in that particular 
  nest() %>% 
  mutate(model01 = purrr::map(data, ~ lm(TotWetHa.nwi ~ TotWetHa.3f, data = .)))
  
# uni.red.ImgYr.nDF <- uni.red.ImgYr.nDF %>% 
#   mutate(model02 = purrr::map(data, ~ lm(TotWetHa.nwi ~ TotWetHa.3f, data = .)))  
# Extract model summaries:
uni.red.ImgYr.nDF <- uni.red.ImgYr.nDF %>% unnest(model01 %>% purrr::map(broom::glance))
uni.red.ImgYr.nDF %>% 
  # ggplot(aes(x = adj.r.squared, y= reorder(IMAGE_YR,adj.r.squared))) +
  # ggplot(aes(x = IMAGE_YR, y = adj.r.squared)) +
  ggplot(aes(x = IMAGE_YR, y = adj.r.squared)) +
  # geom_barh(stat = 'identity', aes(fill = p.value)) +
  geom_bar(stat = 'identity', fill = 'ivory2', color = 'ivory3') +
  geom_point() +
  xlim(1970,2010) +
  labs(y=expression('adjusted r'^{2}), x = "NWI image year", title = "Linear models: 3F area ~ NWI area")

## data table
uni.red.ImgYr.nDF %>% 
  datatable(filter='top')

NA

Climate modeling:

# wetl.uni01red_wide1
vpd.aa <- read_csv("data/Final_shapefile/V15_data/v15centrPRISMNormal_LU.csv")
vpd.aa <- vpd.aa %>% 
  dplyr::select(AAv14id, vpdMx30yr, PrismPPT, PrismTMean) %>% 
  na.omit()
#join with the tabulated full data set
vpd.aa <- left_join(vpd.aa, wetl.uni01red_wide1, by = "AAv14id")
vpd.aa %>% 
  mutate(PercAAWetDiff3f.redNWI = PercWet3f-PercWetNWI) %>% 
  mutate(vpdInt = as.integer(vpdMx30yr)) %>% 
  ggplot(aes(x = vpdMx30yr, y = PercAAWetDiff3f.redNWI)) +
  # geom_point(alpha = 0.05) +
  # geom_density_2d() +
  geom_smooth(method = 'lm') +
  geom_hline(yintercept = 0, color = 'red', alpha=0.41, lty = 'dashed', size = 1.2) +
  scale_y_continuous(limits = quantile(vpd.aa$PercAAWetDiff3f.redNWI, c(0.05, 0.95)))

Unknown column 'PercAAWetDiff3f.redNWI'

### boxplots!
pp1 <- vpd.aa %>% 
  mutate(PercAAWetDiff3f.redNWI = PercWet3f-PercWetNWI) %>% 
  mutate(vpdInt = as.factor(as.integer(vpdMx30yr))) %>% 
  ggplot(aes(x = vpdMx30yr, y = PercAAWetDiff3f.redNWI)) +
  geom_boxplot(aes(group=vpdInt), outlier.shape = NA, fill = "ivory2") +
  labs(x = "Maximum vapor pressure deficit (hPa)", y = "Difference in % AA mapped", title = "VPD") +
  # geom_smooth(method = "lm", se = FALSE) +
  geom_hline(yintercept = 0, color = 'red', alpha=0.41, lty = 'dashed', size = 1.2) +
  theme_smFacet
# ggsave("vpd_AllAA_boxplot.png",width = 6, height = 4.5, dpi = 300)
  
  
## boxplot: ppt
pp2 <- vpd.aa %>% 
  mutate(PercAAWetDiff3f.redNWI = PercWet3f-PercWetNWI) %>% 
  mutate(PPTInt = as.factor(as.integer(PrismPPT))) %>%
  mutate(quartile.PPT = ntile(PrismPPT, 10)) %>% 
  # mutate(t = map_int(.x=.$PrismPPT, .f = cut_equal_ranges))
  # ggplot(aes(x = quartile.PPT, y = PercAAWetDiff3f.redNWI)) +
  ggplot(aes(x = PrismPPT, y = PercAAWetDiff3f.redNWI)) +
  # geom_point(alpha = .051, color = 'darkgrey', size = 2.5) +
  # geom_boxplot(aes(group=t), outlier.shape = NA, fill = "ivory2") +
  labs(x = "Precipitation (mm)", y = "Difference in % AA mapped", title = "Precipitation") +
  geom_smooth(method = "lm", se = FALSE) +
  geom_hline(yintercept = 0, color = 'red', alpha=0.6, lty = 'dashed', size = 1.2) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_smFacet
# ggsave("PrimsPPT_AllAA_boxplot.png",width = 6, height = 4.5, dpi = 300)
## boxplot: temp
pp3 <- vpd.aa %>% 
  mutate(PercAAWetDiff3f.redNWI = PercWet3f-PercWetNWI) %>% 
  mutate(TempInt = as.factor(as.integer(PrismTMean))) %>%
  mutate(quartile.PrismTMean = ntile(PrismTMean, 10)) %>% 
  # mutate(t = map_int(.x=.$PrismPPT, .f = cut_equal_ranges))
  # ggplot(aes(x = quartile.PPT, y = PercAAWetDiff3f.redNWI)) +
  ggplot(aes(x = PrismTMean, y = PercAAWetDiff3f.redNWI)) +
  # geom_point(alpha = .051, color = 'darkgrey', size = 2.5) +
  geom_boxplot(aes(group=TempInt), outlier.shape = NA, fill = "ivory2") +
  labs(x = "Temperature (C)", y = "Difference in % AA mapped", title = "Mean temperature") +
  geom_smooth(method = "lm", se = FALSE) +
  geom_hline(yintercept = 0, color = 'red', alpha=0.6, lty = 'dashed', size = 1.2) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_smFacet
# ggsave("PrimsTEMP_AllAA_boxplot.png",width = 6, height = 4.5, dpi = 300)
# pp3 <- cowplot::plot_grid(pp1, pp2, labels = c("A", "B"), nrow = 2, align = "V")
# save_plot(pp3, filename = "2panel_plot.png",ncol = 1, base_height = 6, base_width = 7.5)
# ### Ternaryplot ... interstingb but not now
# library(ggtern)
# data('Feldspar')
# ggtern(Feldspar,aes(Ab,An,Or)) + 
#  geom_density_tern(aes(color=..level..),bins=5) +
#  geom_point()
# 
# ggtern(vpd.aa,aes(PrismPPT,PrismTMean,vpdMx30yr)) + 
#  geom_density_tern(aes(color=..level..),bins=5) +
#  geom_point()
# 
# # ### hex plot
# # vpd.aa %>% 
# #   mutate(PercAAWetDiff3f.redNWI = PercWet3f-PercWetNWI) %>% 
# #   mutate(vpdInt = as.factor(as.integer(vpdMx30yr))) %>% 
# #   ggplot(aes(x = vpdMx30yr, y = PercAAWetDiff3f.redNWI)) +
# #   geom_hex(bins=35, alpha = 0.95) +
# #   scale_fill_viridis(direction = -1) +
# #   geom_quantile() 

Climatevars scatterplots

plot01 <- vpd.aa %>%
  ggplot(aes(x = PrismPPT, y = diff.AApercWet)) +
  geom_point(alpha = 0.08) +
  geom_smooth() +
  geom_hline(yintercept = 0, lty = 'dashed',color='red',size=1.2) +
  labs(x="Precipitation (mm)", y = "Difference in %AA mapped")
plot01

# save_plot(plot = plot01,filename = "ppt_vs_DiffAA.png")

Looking at 3f polygon ditribution

F3poly <- map2(.x = "data/Final_shapefile/v15_data", .y = "Fin3F_AAv14DIS", .f =  readOGR) ### Note: this also works

OGR data source with driver: ESRI Shapefile 
Source: "data/Final_shapefile/v15_data", layer: "Fin3F_AAv14DIS"
with 29490 features
It has 6 fields

Z-dimension discarded

# ### View attribute data
F3poly.tbl <- as_tibble(F3poly[[1]]@data) 
# names(F3poly.tbl)
F3poly.tbl <- left_join(F3poly.tbl,ecoreg.lu,by="AAv14id")

joining character vector and factor, coercing into character vector

# F3poly.tbl %>% datatable()
## summarize total AA area
F3poly.tbl %>% 
  # names() %>% 
  summarize(sum.polyHa = sum(PolyAreaHa), mean.polyHa = mean(PolyAreaHa, na.rm = TRUE), med.polyHa = median(PolyAreaHa), cnt.poly = n(), min.polyHa = min(PolyAreaHa), max.polyHa = max(PolyAreaHa)) %>%
  mutate(sum.polyAcres = 2.47105*sum.polyHa) %>% 
  mutate(mean.polyAcres = 2.47105*mean.polyHa) %>% 
  mutate(med.polyAcres = 2.47105*med.polyHa) %>% 
  mutate(min.polyAcres = 2.47105*min.polyHa) %>% 
  mutate(max.polyAcres = 2.47105*max.polyHa)

F3poly.tbl <- F3poly.tbl %>%
  mutate(PolyAreaAcre = PolyAreaHa*2.47105) 
  
F3poly.tbl %>% datatable()

It seems your data is too big for client-side DataTables. You may consider server-side processing: http://rstudio.github.io/DT/server.htmlIt seems your data is too big for client-side DataTables. You may consider server-side processing: http://rstudio.github.io/DT/server.html