# #Here's a page about accessing this data via Socrata: https://dev.socrata.com/foundry/opendata.maryland.gov/ed4q-f8tm
# #Here's a guide to clauses for querying: https://dev.socrata.com/docs/queries/
# 
# #For this example, I am looking for counties that have aquaculture along Maryland's Chesapeake Bay. I would like to know the assessed value of waterfront (PFLW or Field #65 = 1 (waterfront)) parcels with structures on them (NFMIMPVL - improved value; NFMTTLVL - total value of parcel land + structure), and what the pacel is zoned for (DESCLU - land use). 
# 
# ####Background Connect to server####
# water.view.parcels <- read.socrata("https://opendata.maryland.gov/resource/ed4q-f8tm.json?property_factors_location_waterfront_mdp_field_pflw_sdat_field_65=Water View (2)")
# 
# water.parcels <- read.socrata("https://opendata.maryland.gov/resource/ed4q-f8tm.json?$where=starts_with(property_factors_location_waterfront_mdp_field_pflw_sdat_field_65, 'Water')")
# colnames(water.parcels)
# 
# #subset my columns
# needed.cols <- water.parcels %>% select(jurisdiction_code_mdp_field_jurscode, county_name_mdp_field_cntyname, mdp_longitude_mdp_field_digxcord_converted_to_wgs84, mdp_latitude_mdp_field_digycord_converted_to_wgs84, land_use_code_mdp_field_lu_desclu_sdat_field_50, property_factors_location_waterfront_mdp_field_pflw_sdat_field_65, base_cycle_data_date_assessed_yyyy_mm_sdat_field_158, current_cycle_data_improvements_value_mdp_field_names_nfmimpvl_curimpvl_and_salimpvl_sdat_field_165, current_cycle_data_land_value_mdp_field_names_nfmlndvl_curlndvl_and_sallndvl_sdat_field_164, prior_assessment_year_total_assessment_sdat_field_161, current_assessment_year_total_assessment_sdat_field_172)
#   
# #select bay counties  
# bay.parcels <- subset(needed.cols, jurisdiction_code_mdp_field_jurscode %in% c("ANNE","CALV", "STMA", "CHAR", "KENT", "QUEE", "TALB", "DORC", "SOME", "WICO"))
# 
# #rename columns
# bay.data <- bay.parcels %>% rename(
#   county.code  = jurisdiction_code_mdp_field_jurscode,
#   county.name = county_name_mdp_field_cntyname,
#   long = mdp_longitude_mdp_field_digxcord_converted_to_wgs84,
#   lat = mdp_latitude_mdp_field_digycord_converted_to_wgs84,
#   land.use = land_use_code_mdp_field_lu_desclu_sdat_field_50,
#   waterfront = property_factors_location_waterfront_mdp_field_pflw_sdat_field_65,
#   date.assessed = base_cycle_data_date_assessed_yyyy_mm_sdat_field_158,
#   improvement.value = current_cycle_data_improvements_value_mdp_field_names_nfmimpvl_curimpvl_and_salimpvl_sdat_field_165,
#   land.value = current_cycle_data_land_value_mdp_field_names_nfmlndvl_curlndvl_and_sallndvl_sdat_field_164,
#   prior_assessment_year_total_assessment = prior_assessment_year_total_assessment_sdat_field_161,
#   current_assessment_year_total_assessment =current_assessment_year_total_assessment_sdat_field_172)
# bay.data <- data.frame(bay.data)
# 
# #Create a Shapefile from the data
# bay.data$lat <- as.numeric(bay.data$lat)
# bay.data$long <-  as.numeric(bay.data$long)
# bay.data <- bay.data[-which(is.na(bay.data$lat)),]
# coordinates(bay.data) <- c("long", "lat") 
# crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")  
# crs(bay.data) <- crs.geo  
# plot(bay.data)
# 
# ####From Here: Export the untransformed data as a shapefile####
# outfile <- "bay.houses.shp"
# shapefile(bay.data, outfile, overwrite=TRUE)
bay.data <- shapefile("bay.houses.shp")
bay.data <- spTransform(bay.data, CRSobj =  "+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs")
bay.data <- st_as_sf(bay.data)


#Bring in shoreline and buffer it
setwd("~/Desktop/MEES Graduate School Work/Thesis Work/Data")
shoreline <- shapefile("shoreline_line.shp")
shoreline <- spTransform(shoreline, CRSobj =  "+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs")
shoreline <- st_as_sf(shoreline)
shoreline_simp = st_simplify(shoreline, dTolerance = 1000)  # 1000 m
#plot(shoreline_simp)
shoreline_buff_50m = st_buffer(shoreline_simp, dist = 50)
#plot(shoreline_buff_50m)
#houses in 50m of shoreline
houses <- st_intersection(bay.data, shoreline_buff_50m)
drops <- c("imprvm_","land_vl", "prr____") # list of col names
houses <- houses[,!(names(houses) %in% drops)] #remove columns 
houses.xy <- st_coordinates(houses)
houses.df <- as.data.frame(houses)
CPI <- read.csv("Average Annual CPI.csv")
cpi <- as.data.frame(CPI)
houses.df$dt_ssss <- substr(houses.df$dt_ssss, 0, 4)
houses.df$dt_ssss <- as.numeric(houses.df$dt_ssss)
houses.df$crr____ <- as.numeric(houses.df$crr____)


head(cpi)
year.rate <- cpi[, c("Year", "multiplier_2020")]
colnames(year.rate)[1] <- "dt_ssss"
year.rate$dt_ssss <- as.numeric(year.rate$dt_ssss)

adjusted.houses.df <- houses.df %>%
  left_join(year.rate, by = c("dt_ssss"="dt_ssss"))%>%
  mutate(value_2020 = crr____* multiplier_2020)
adjusted.houses.df <- adjusted.houses.df[-c(7:9, 11)]
adjusted.houses.df <- st_as_sf(adjusted.houses.df)
adjusted.houses.df <- as(adjusted.houses.df, "Spatial")
writeOGR(adjusted.houses.df, ".", "adjusted.houses.df", driver = "ESRI Shapefile", overwrite_layer = T)

#buffer houses by significant distance values found in Stump thesis
adjusted.houses.df <- st_as_sf(adjusted.houses.df)
houses_buff_100m = st_buffer(adjusted.houses.df, dist = 100)
houses_buff_100m <- houses_buff_100m %>%
  add_column(Distance = "100")

houses_buff_300m = st_buffer(adjusted.houses.df, dist = 300)
houses_buff_300m <- houses_buff_300m %>%
  add_column(Distance = "300")

houses_buff_400m = st_buffer(adjusted.houses.df, dist = 400)
houses_buff_400m <- houses_buff_400m %>%
  add_column(Distance = "400")

houses_buff_500m = st_buffer(adjusted.houses.df, dist = 500)
houses_buff_500m <- houses_buff_500m %>%
  add_column(Distance = "500")


# houses_buff_500m <- as(houses_buff_500m, "Spatial")
# writeOGR(houses_buff_500m, ".", "houses_buff_500m", driver = "ESRI Shapefile", overwrite_layer = T)
# 
# houses_buff_400m <- as(houses_buff_400m, "Spatial")
# writeOGR(houses_buff_400m, ".", "houses_buff_400m", driver = "ESRI Shapefile", overwrite_layer = T)
# 
# houses_buff_300m <- as(houses_buff_300m, "Spatial")
# writeOGR(houses_buff_300m, ".", "houses_buff_300m", driver = "ESRI Shapefile", overwrite_layer = T)
# 
# houses_buff_100m <- as(houses_buff_100m, "Spatial")
# writeOGR(houses_buff_100m, ".", "houses_buff_100m", driver = "ESRI Shapefile", overwrite_layer = T)

####background information for looping####
#avg lease size is 16.56 acres regardless of gear type // 19.71 acres for submerged // 5.6 for wc
#1. calculate water acres of buffer per parcel
#2. randomly place leases
#3. calculate % of area of each buffer per parcel consumed by leases
#4. calculate chanage in parcel value based on lease type, space consumed, and distance 

cb <- shapefile("Chesapeake_Bay_92_Segments.shp")
cb <- spTransform(cb, CRSobj =  "+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs")

cb_sf <- st_as_sf(cb)
cb_simp = st_simplify(cb_sf, dTolerance = 50)  # 50 m
bay_side_100m <- st_intersection(cb_simp, houses_buff_100m)
bay_side_300m <- st_intersection(cb_simp, houses_buff_300m)
bay_side_400m <- st_intersection(cb_simp, houses_buff_400m)
bay_side_500m <- st_intersection(cb_simp, houses_buff_500m)


plot(st_geometry(bay_side_100m[10]), col = sf.colors(12, categorical = TRUE), border = 'grey',
     axes = F)
plot(st_geometry(bay_side_300m[10]), col = sf.colors(10, categorical = TRUE), add=T)
plot(st_geometry(bay_side_400m[10]), col = sf.colors(8, categorical = TRUE), add=T)
plot(st_geometry(bay_side_500m[10]), col = sf.colors(6, categorical = TRUE), add=T)


####extra treatments####
#Log transform value data
bay.data$improvement.value <- as.numeric(bay.data$improvement.value)
bay.data$improvement.value <- log(bay.data$improvement.value)

bay.data$land.value <- as.numeric(bay.data$land.value)
bay.data$land.value <- log(bay.data$land.value)
is.na(bay.data) <- sapply(bay.data, is.infinite)
bay.data <- na.omit(bay.data)

#Print a histogram by county
v <- count(bay.data, county.name,improvement.value)
print_list <- unique(bay.data$county.name)
for (i in seq_along(print_list)) {
  df=subset(v, county.name==print_list[i])
  plt <- ggplot(df, aes(improvement.value, n)) +
    geom_bar(stat = "identity", width = 0.9) +
    labs(y = "Property Value", x = "") + 
    xlim(0,25) +
    ylim(0,200)+
    labs(caption = print_list[i]) +
    theme_void()
  b <- plt + theme(plot.caption = element_text(hjust=0.5, face="bold", size = 20))
  
  png(paste0(print_list[i], ".png"), 4, 4, "in", res = 100, bg="transparent")
  print(b)
  dev.off()
}


####Houses with 500m buffer, look at how many leases already impact house values####
adjusted.houses.df <- shapefile("Houses/adjusted.houses.df.shp")
adjusted.houses.df <- spTransform(adjusted.houses.df, CRSobj =  "+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs")
adjusted.houses.df <- st_as_sf(adjusted.houses.df)
adjusted.houses.df <- subset(adjusted.houses.df, land_us=="Residential (R)")

#residential houses with 500m buffer
houses_sf_500m <- st_buffer(adjusted.houses.df, 500)

#residential houses with unique row numbers
adjusted.houses.df <- adjusted.houses.df%>% dplyr::mutate(ID = row_number())
#remove geometry 
adjusted.houses.df <- st_set_geometry(adjusted.houses.df, NULL)
adjusted.houses.df$cnty_nm <- as.factor(adjusted.houses.df$cnty_nm)
count_table <- adjusted.houses.df %>%
  group_by(cnty_nm)%>%
  count()
View(count_table)
sum(adjusted.houses.df$value_2020,na.rm=TRUE)

leases <- shapefile("leases/oysterleases.shp")
leases <- spTransform(leases, CRSobj =  "+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs")
leases.sf <- st_as_sf(leases)

houses_sf_500m <- st_intersection(leases.sf, houses_sf_500m)

houses_sf_21leases <- st_join(house_value_sf,houses_sf_500m)
houses_sf_21leases <- houses_sf_21leases%>%na.omit(L0Active_O)


#houses_sf_500m <- st_set_geometry(houses_sf_500m, NULL)
#houses_sf_500m <- as.data.frame(houses_sf_500m)




####Calculating policy impacts ####
#NEED TO CHECK IF THERE IS A LEASE THAT OVERLAPS AT MULTIPLE DISTANCES, SELECT MAX DISTANCE FOR CALCULATIONS

house_value_impacts <- shapefile("house_values/a_prcnt_area_bygearscenario.shp")
house_value_sf <- st_as_sf(house_value_impacts)

multipliers <- read.csv("Stump_values.csv")

#Rename column in dataframe
names(multipliers)[1] <- "BUFF_DIST"

#Rename categorical value in dataframe
multipliers$gear <- if_else(multipliers$gear == "sub",
         "submerged",
         "wc")
#join value table based on conditions
house_value_sf <- left_join(house_value_sf, multipliers, by = c('BUFF_DIST', 'gear')) 
house_value_sf$value_2020 <- as.numeric(house_value_sf$value_2020)
house_value_sf$mult <- as.numeric(house_value_sf$mult)
house_value_sf$prcnt_over <- as.numeric(house_value_sf$prcnt_over)
house_impacted_value_cal_nogeom <- st_set_geometry(house_value_sf, NULL)

#bring in nicknames CSV
nicknames <- read.csv("policy name book.csv")
house_impacted_value_cal_nogeom <- left_join(house_impacted_value_cal_nogeom, nicknames, by = 'scenari')
house_impacted_value_cal_nogeom <- house_impacted_value_cal_nogeom %>%
  na.omit(mult)

View(house_impacted_value_cal_nogeom)
####Calculations####

#calculate value change (pcnt_change is percent overlap * change value, policy value is new value of parcel, net_value_change is $ value removed from parcel)
house_impacted_value_cal <- house_impacted_value_cal_nogeom %>%
  drop_na(gear)%>%
  mutate(pcnt_change=prcnt_over*mult, 
         net_value=value_2020+((pcnt_change/100)*value_2020), 
         net_value_change= ((pcnt_change/100)*value_2020)
         # value_change_Policy = ifelse(nickname !='Baseline Conditions (2018)',
         # net_value-net_value[nickname == 'Baseline Conditions (2018)'],
         #                              0),
         # policy_value = value_2020 + value_change_Policy,
         #policy_pcnt_change = (policy_value/value_2020)*100
         ) 
View(house_impacted_value_cal)
#Now, for each OBJECTID and nickname, select the max pcnt_change
house_impacted_value_cal <- house_impacted_value_cal %>%
  group_by(OBJECTID, nickname)%>%
  filter(pcnt_change == max(pcnt_change, na.rm=TRUE))

#% value changes by policy
#net_value_change is ∆ in $ from 2020 to 2036, net value is the 2036 value, sum_net_value_change is the average ∆$s for properties change when 2018 conditions are controlled for, sum_policy_value is the value of properties when 2018 conditions are controlled for, mean_policy_pcnt_change = % house value ∆s because of policy with 2018 controlled for
county_impacted_value_calc <- house_impacted_value_cal %>%
  drop_na(gear) %>%
  group_by(cnty_nm, nickname) %>%
  summarise(mean_overlapping_area = mean(prcnt_over),
            sum_value_2020 = sum(value_2020),
            sum_net_value_change = sum(net_value_change))%>%
#            sum_net_value_change_change = sum(net_value_change),
#            sum_net_value_change = sum(value_change_Policy),
#            sum_policy_value = sum(policy_value), 
#            mean_policy_pcnt_change = mean(policy_pcnt_change)
  mutate(percent_delta = (sum_net_value_change/sum_value_2020)*100
    )
View(county_impacted_value_calc)

plot_valuechange <- county_impacted_value_calc %>%
  drop_na(cnty_nm)%>%
#  drop_na(mean_policy_pcnt_change)%>%
  group_by(nickname)%>%
  ggplot(aes(x=reorder_within(cnty_nm, percent_delta, nickname), 
             y=percent_delta, fill = cnty_nm)) +
  geom_bar(stat='identity') +
  labs(fill = "County") +
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())+
  facet_wrap(~reorder(nickname, percent_delta), scales = "free_x")+
  ylab("Total Change in R Waterfront Home Value (%)")+
  scale_y_continuous(limits=c(-30, NA))
plot(plot_valuechange)

# $$$value changes by policy
plot_valuechange <- county_impacted_value_calc %>%
  drop_na(cnty_nm)%>%
  group_by(nickname)%>%
  mutate(group_mean = mean(sum_net_value_change)) %>%
  ggplot(aes(x=reorder_within(cnty_nm, sum_net_value_change, nickname), 
             y=sum_net_value_change, fill = cnty_nm)) +
  geom_bar(stat='identity') +
  labs(fill = "County") +
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())+
  facet_wrap(~reorder(nickname, group_mean), scales = "free_x")+
  ylab("Total Change in R Waterfront Home Value (2020 US $)")+
  scale_y_continuous(labels = comma)+
  geom_hline(aes(yintercept=group_mean), colour="grey70")
plot(plot_valuechange)

write.csv(plot_valuechange,"Impact of policy on house values.csv", row.names = FALSE)


#statewide average by scenario (change to look at median vs mean)
state_valuechange <- county_impacted_value_calc %>%
  group_by(nickname) %>%
  summarise(pcnt_val_change=sum(sum_net_value_change)/sum(sum_value_2020))
View(state_valuechange)

#Which counties have the greatest % houses impacted per scenario? count of houses?

house_count <- house_impacted_value_cal_nogeom%>% 
  group_by(cnty_nm, nickname, BUFF_DIST)%>%
  summarise(count=n_distinct(OBJECTID))
house_count$count <- as.numeric(house_count$count)
View(house_count)
house_count <- left_join(count_table, house_count, by = c('cnty_nm'))

house_count_by_buffer <- house_count%>%
  group_by(cnty_nm, nickname)%>%
  mutate(perc= ((count/n)*100)) %>%
  na.omit(nickname)
View(house_count)

#House count regardless of buffer at which lease falls
house_count <- house_impacted_value_cal_nogeom%>% 
  group_by(cnty_nm, nickname)%>%
  summarise(count=n_distinct(OBJECTID))
house_count$count <- as.numeric(house_count$count)
house_count <- left_join(count_table, house_count, by = c('cnty_nm'))

house_count <- house_count%>%
  group_by(cnty_nm, nickname)%>%
  mutate(perc= ((count/n)*100)) %>%
  na.omit(nickname)
View(house_count)

house_count <- house_count%>%
  group_by(cnty_nm, nickname)%>%
  arrange(desc(perc), .by_group = TRUE)
write.csv(house_count,"Count of houses impacted by policies.csv", row.names = FALSE)

plot_percHousebyScenario <- house_count %>%
  ggplot(aes(x=reorder_within(cnty_nm, -perc, nickname), y=perc, fill = cnty_nm)) +
  geom_bar(stat='identity') +
  labs(fill = "County") +
  geom_text(aes(label=count), vjust=-0.3, size=2.5)+
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())+
  ylab("Percent of Waterfront Residential Houses Impacted") +
  ylim(0,100)+
  facet_wrap(~reorder(nickname, -perc), scales = "free_x")
plot_percHousebyScenario