From 9e8f600ddfb3ed40fdbe2fb137324d05dafd9e58 Mon Sep 17 00:00:00 2001 From: Alex Date: Fri, 6 Mar 2026 14:31:39 -0700 Subject: [PATCH] Intial cleanup of git data --- .gitignore | 2 + 1_Data_Proc.r | 258 +++++++++++++++++++++++++++++++++++ Hydrobase.r => 2_Hydrobase.r | 0 Data/Readme_Data.md | 5 + Scripts/Distance_Functions.r | 75 ++++++++++ 5 files changed, 340 insertions(+) create mode 100755 1_Data_Proc.r rename Hydrobase.r => 2_Hydrobase.r (100%) mode change 100644 => 100755 create mode 100644 Data/Readme_Data.md create mode 100644 Scripts/Distance_Functions.r diff --git a/.gitignore b/.gitignore index f26b794..9d14b0c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ # ---> R +Data/Input_Data/* +Data/Input_Data.zip # History files .Rhistory .Rapp.history diff --git a/1_Data_Proc.r b/1_Data_Proc.r new file mode 100755 index 0000000..701a93d --- /dev/null +++ b/1_Data_Proc.r @@ -0,0 +1,258 @@ +library(tidyverse) +library(fixest) +#function "i" needed for a factor +crop <- read_csv("./Data/Input_Data/CROPS.csv",col_types = cols(PARCEL_ID=col_character(),SW_WDID1=col_character(),SW_WDID2=col_character(),SW_WDID3=col_character(),SW_WDID4=col_character(),SW_WDID5=col_character(),GW_ID1=col_character(),GW_ID2=col_character(),GW_ID3=col_character(),GW_ID4=col_character(),GW_ID5=col_character(),GW_ID6=col_character(),GW_ID7=col_character(),GW_ID6=col_character(),GW_CLASS6=col_character(),GW_ID7=col_character(),GW_CLASS7=col_character(),GW_CLASS1=col_character(),GW_CLASS2 =col_character(),GW_CLASS3 =col_character(),GW_CLASS4 =col_character(),GW_CLASS5 =col_character(),GW_CLASS6 =col_character(),GW_CLASS7 =col_character() )) %>% select(PARCEL_ID,Year=CAL_YEAR,CROP_TYPE,IRR_TY=IRRIG_TYPE,CROP_AREA,GW_ID1,GW_ID2,GW_ID3,GW_ID4,GW_ID5,GW_ID6,GW_ID7,SW_WDID1,SW_WDID2,SW_WDID3,SW_WDID4,SW_WDID5) + +crop$IRR_TY <- ifelse(crop$IRR_TY=='FLOOD','F','S') + +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="NEW_ALFALFA","ALFALFA",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="BARLEY","SMALL_GRAINS",crop$CROP_TYPE) +#blue grass, cover crop, sorghum grain, vegetables, wheat, fall wheat, spring wheat, and corn. +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="COVER_CROP","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="VEGETABLES","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="SORGHUM_GRAIN","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="WHEAT","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="WHEAT_FALL","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="WHEAT_SPRING","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="CORN","OTHER",crop$CROP_TYPE) +crop$CROP_TYPE <- ifelse(crop$CROP_TYPE=="","OTHER",crop$CROP_TYPE) +crop$POT <- ifelse(crop$CROP_TYPE=="POTATOES",crop$CROP_AREA,0) +crop$SG<- ifelse(crop$CROP_TYPE=="SMALL_GRAINS",crop$CROP_AREA,0) +crop$ALF <- ifelse(crop$CROP_TYPE=="ALFALFA",crop$CROP_AREA,0) +crop$OT <- ifelse(crop$CROP_TYPE=="OTHER",crop$CROP_AREA,0) +crop$GP <- ifelse(crop$CROP_TYPE=="GRASS_PASTURE",crop$CROP_AREA,0) + +SW_LIST <- crop %>% group_by(SW_WDID1) %>% filter(!is.na(SW_WDID1)) %>% summarize(num=n()) %>% arrange(desc(num)) %>% slice(1:10) %>% pull(SW_WDID1) +SW_COL <- paste0("SW_",SW_LIST) +SW_DUM <- function(X){ +return(ifelse(crop$SW_WDID1 %in% X | crop$SW_WDID2 %in% X | crop$SW_WDID3 %in% X | crop$SW_WDID3 %in% X | crop$SW_WDID5 %in% X,1,0)) +} +crop$SW_2000812 <- SW_DUM(SW_LIST[1]) +crop$SW_2000623 <- SW_DUM(SW_LIST[2]) +crop$SW_2000631 <- SW_DUM(SW_LIST[3]) +crop$SW_2000753 <- SW_DUM(SW_LIST[4]) +crop$SW_2200593 <- SW_DUM(SW_LIST[5]) +crop$SW_3500570 <- SW_DUM(SW_LIST[6]) +crop$SW_2000829 <- SW_DUM(SW_LIST[7]) +crop$SW_2000798 <- SW_DUM(SW_LIST[8]) +crop$SW_2100521 <- SW_DUM(SW_LIST[9]) +crop$SW_2200541 <- SW_DUM(SW_LIST[10]) + + + +wells <- read_rds("./Data/Input_Data/Well_Pumping.rds") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID1,all_of(SW_COL)) %>% distinct +m1 <- wells %>% inner_join(mt,multiple="all") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID2,all_of(SW_COL)) %>% distinct +m2 <- wells %>% inner_join(mt,multiple="all") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID3,all_of(SW_COL)) %>% distinct +m3 <- wells %>% inner_join(mt,multiple="all") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID4,all_of(SW_COL)) %>% distinct +m4 <- wells %>% inner_join(mt,multiple="all") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID5,all_of(SW_COL)) %>% distinct +m5 <- wells %>% inner_join(mt,multiple="all") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID6,all_of(SW_COL)) %>% distinct +m6 <- wells %>% inner_join(mt,multiple="all") +mt <- crop %>% select(PARCEL_ID,Year=Year,wdid=GW_ID7,all_of(SW_COL)) %>% distinct +m7 <- wells %>% inner_join(mt,multiple="all") +m <- rbind(m1,m2,m3,m4,m5,m6,m7) %>% rename(GW_wdid=wdid) + +crop <- crop %>% select(PARCEL_ID,Year,CROP_TYPE,CROP_AREA,IRR_TY,POT,SG,ALF,GP,OT) +#Crep + +CREP_PERM_2014 <- c(2013956,2005950,2005951,2005955,2014091,2006322,2006321,2014092, 2006327,2006328,2006325,2006326,2006332,2006331,2006684,2006685,2006686,2014274,2014107,2005512,2005448) +CREP_PERM_2015 <- c(2005098,2008177,2008178,2013955,2006005,2006656,2005171,2006655,2006323,2006324,2014088,2008223,2008224,2008225,2014054,2705248) +CREP_PERM_2016 <- c(2005476,2005537,2005538,2014266,2005769,2005770,2005771,2014270,2005766,2005767,2005768,2014267,2014268,2008439,2008440,2208441,2009992,2009197,2014045,2014046,2006003,2006004,2006653,2006654,2014311,2005952,2005953,2005954,2005121,2008772) +CREP_PERM_2018 <- c(2014309,2006334,2014080,2006333,2006337,2006338,2014081,2014082) +CREP_PERM_2020 <- c(2010499,2010500,2013906,2014178,2006528,2006529,2014189,2009110,2008239,2705472) +CREP_PERM_2021 <- c(2008146,2008147,2013953,2005665,2005664,2005666,2005662,2005663,2010786,2010787,2008389,2014023) +CREP_PERM_2014 <- cbind(CREP_PERM_2014,rep(2014,length(CREP_PERM_2014)),rep(1,length(CREP_PERM_2014))) +CREP_PERM_2015 <- cbind(CREP_PERM_2015,rep(2015,length(CREP_PERM_2015)),rep(1,length(CREP_PERM_2015))) +CREP_PERM_2016 <- cbind(CREP_PERM_2016,rep(2016,length(CREP_PERM_2016)),rep(1,length(CREP_PERM_2016))) +CREP_PERM_2018 <- cbind(CREP_PERM_2018,rep(2018,length(CREP_PERM_2018)),rep(1,length(CREP_PERM_2018))) +CREP_PERM_2020 <- cbind(CREP_PERM_2020,rep(2020,length(CREP_PERM_2020)),rep(1,length(CREP_PERM_2020))) +CREP_PERM_2021 <- cbind(CREP_PERM_2021,rep(2021,length(CREP_PERM_2021)),rep(1,length(CREP_PERM_2021))) + +CREP_PERM <- rbind(CREP_PERM_2014,CREP_PERM_2015,CREP_PERM_2016,CREP_PERM_2018,CREP_PERM_2020,CREP_PERM_2021) +colnames(CREP_PERM) <- c("GW_wdid","CREP_Year","PERM_CREP") + + + +CREP_TEMP_2014 <- c(2005642,2005643,2014474,2006153,2013962,2006478,2008677,2008678,2012887,2005857,2008391,2705126,2705519,2706148,2705342,2706196,2705341,2706195) +CREP_TEMP_2015 <- c(2008155,2008156,2008129,2008130,2014244,2005921,2005941,2006283,2006525,2006335,2006336,2014086,2014087,2705318,2012537,2014288,2706253,2705317,2705186,2705328,2705054) +CREP_TEMP_2016 <- c(2005133,2005533,2005886,2005868,2005595,2013377,2013618,2706014,2706246) +CREP_TEMP_2017 <- c(2009113,2705344,2705067,2705068,2705523,2705069,2705070) +CREP_TEMP_2018 <- c(2005774,2005775,2705293,2705290,2705224,2705225,2705225,2705224,2705197,2705359,2705246,2706237,2705006,2705790,2705184,2705185,2705356,2705327) +CREP_TEMP_2019 <- c(2006376,2006375,2005127,2005168,2013784,2705259,2705021,2705020,2706194) +CREP_TEMP_2020 <- c(2705347,2706258) +CREP_TEMP_2022 <- c(2006678,2006679,2005923,2005116) +CREP_TEMP_2023 <- c(2009084,2009082,2013887,2009940,2013577,2005038,2010590,2014490,2705247,2705080) + +CREP_TEMP_2014 <- cbind(CREP_TEMP_2014,rep(2014,length(CREP_TEMP_2014)),rep(0,length(CREP_TEMP_2014))) +CREP_TEMP_2015 <- cbind(CREP_TEMP_2015,rep(2015,length(CREP_TEMP_2015)),rep(0,length(CREP_TEMP_2015))) +CREP_TEMP_2016 <- cbind(CREP_TEMP_2016,rep(2016,length(CREP_TEMP_2016)),rep(0,length(CREP_TEMP_2016))) +CREP_TEMP_2017 <- cbind(CREP_TEMP_2017,rep(2017,length(CREP_TEMP_2017)),rep(0,length(CREP_TEMP_2017))) +CREP_TEMP_2018 <- cbind(CREP_TEMP_2018,rep(2018,length(CREP_TEMP_2018)),rep(0,length(CREP_TEMP_2018))) +CREP_TEMP_2019 <- cbind(CREP_TEMP_2019,rep(2019,length(CREP_TEMP_2019)),rep(0,length(CREP_TEMP_2019))) +CREP_TEMP_2020 <- cbind(CREP_TEMP_2020,rep(2020,length(CREP_TEMP_2020)),rep(0,length(CREP_TEMP_2020))) +#CREP_TEMP_2022 <- cbind(CREP_TEMP_2022,rep(2022,length(CREP_TEMP_2022)),rep(0,length(CREP_TEMP_2022))) +#CREP_TEMP_2023 <- cbind(CREP_TEMP_2023,rep(2023,length(CREP_TEMP_2023)),rep(0,length(CREP_TEMP_2023))) +CREP_TEMP <- rbind(CREP_TEMP_2014,CREP_TEMP_2015,CREP_TEMP_2016,CREP_TEMP_2017,CREP_TEMP_2018,CREP_TEMP_2019,CREP_TEMP_2020) +colnames(CREP_TEMP) <- c("GW_wdid","CREP_Year","PERM_CREP") +IN_CREP <- rbind(CREP_TEMP,CREP_PERM) %>% as_tibble %>% mutate(GW_wdid=as.character(GW_wdid)) + +FALL_2020 <- c(2013693,2005731,2013884,2706159,2705471,2705395,2705517) +FALL_2020 <- cbind(FALL_2020,rep(2020,length(FALL_2020)),rep(1,length(FALL_2020))) +FALL_2021 <- c(2008226,2014273,2008591,2009461,2009462,2009114,2010568,2005035,2014173) +FALL_2021 <- cbind(FALL_2021,rep(2021,length(FALL_2021)),rep(1,length(FALL_2021))) +FALL_PROG <- rbind(FALL_2020,FALL_2021) %>% as_tibble +colnames(FALL_PROG) <- c("GW_wdid","FALL_PROG_YEAR","FALL_PROG") + +#Legal Parcel +LP <- read_csv("./Data/Input_Data/LEGAL.csv") %>% rename(County=Layer) %>% select(-fid) +LP$County <- ifelse(is.na(LP$County),"Saguache",LP$County) +LP$County<- ifelse(LP$County=="Saguache",'SG',LP$County) +LP$County<- ifelse(LP$County=="Saguache Parcels",'SG',LP$County) +LP$County<- ifelse(LP$County=="Alamosa Parcels",'AL',LP$County) +LP$County <- ifelse(LP$County=="Rio Grande Parcels",'RG',LP$County) +LP$County <- ifelse(LP$County=="Conejos Parcels",'CN',LP$County) +#Crop to Legal Parcel +CL <- read_csv("./Data/Input_Data/CROP_LEGAL.csv",col_types=cols(PARCEL_ID=col_character())) %>% select(Name,PARCEL_ID,Year=CAL_YEAR) + + +SL <- read_csv("Data/Input_Data/SBD_LEGAL.csv") +SL$SBD1 <- ifelse(SL$SBD==1,1,0) +SL$SBD2 <- ifelse(SL$SBD==2,1,0) +SL$SBD3 <- ifelse(SL$SBD==3,1,0) +SL$SBD4 <- ifelse(SL$SBD==4,1,0) +SL$SBD5 <- ifelse(SL$SBD==5,1,0) +SL$SBD6 <- ifelse(SL$SBD==6,1,0) +SL <- SL %>% select(-SBD) %>% group_by(Name) %>% mutate(SBD1=max(SBD1),SBD2=max(SBD2),SBD3=max(SBD3),SBD4=max(SBD4),SBD5=max(SBD5),SBD6=max(SBD6)) %>% distinct +SL <- LP %>% select(Name) %>% distinct %>% left_join(SL) %>% mutate(across(everything(), ~replace_na(.x, 0))) +LP <- LP %>% left_join(SL) + +CS <- read_csv("Data/Input_Data/SBD_CROP.csv",col_types=cols(PARCEL_ID=col_character())) %>% rename(Year=CAL_YEAR) %>% distinct +CS$SBD1 <- ifelse(CS$SBD==1,1,0) +CS$SBD2 <- ifelse(CS$SBD==2,1,0) +CS$SBD3 <- ifelse(CS$SBD==3,1,0) +CS$SBD4 <- ifelse(CS$SBD==4,1,0) +CS$SBD5 <- ifelse(CS$SBD==5,1,0) +CS$SBD6 <- ifelse(CS$SBD==6,1,0) +CS <- CS %>% select(-SBD) %>% group_by(PARCEL_ID,Year) %>% mutate(SBD1=max(SBD1),SBD2=max(SBD2),SBD3=max(SBD3),SBD4=max(SBD4),SBD5=max(SBD5),SBD6=max(SBD6)) %>% distinct +CS <- crop %>% select(PARCEL_ID,Year) %>% left_join(CS) %>% mutate(across(everything(), ~replace_na(.x, 0))) +m <- m %>% group_by(GW_wdid) %>% mutate(SW_2000812=max(SW_2000812),SW_2000623=max(SW_2000623),SW_2000631=max(SW_2000631),SW_2000753=max(SW_2000753),SW_2200593=max(SW_2200593),SW_3500570=max(SW_3500570),SW_2000829=max(SW_2000829),SW_2000798=max(SW_2000798),SW_2100521=max(SW_2100521),SW_2200541=max(SW_2200541)) %>% ungroup +WSW <- m %>% select(-Year,-AF,-PARCEL_ID) %>% unique +GW_CREP <- m %>% left_join(CL,multiple="all") %>% left_join(IN_CREP) +GW_CREP$IN_CREP <- ifelse(is.na(GW_CREP$CREP_Year),0,1) + +#Well data set +well_df <- GW_CREP %>% left_join(m) +well_df <- well_df %>% full_join(crop) +well_df <- well_df %>% left_join(CS) +#well_df %>% filter(SBD1==0,SBD2==0,SBD3==0,SBD4==0,SBD5==0,SBD6==0) +well_df[which(is.na(well_df$GW_wdid)),1] <- "NONE" +well_df$AF<- ifelse(is.na(well_df$AF),0,well_df$AF) +#One point is -0.125 checked Hydrobase to confirm values. Either bad entry or augmentation +well_df[well_df$AF<0,]$AF <- 0 +#well_df$IN_CREP<- ifelse(is.na(well_df$IN_CREP),0,well_df$IN_CREP) +well_df$IN_CREP <- well_df$GW_wdid %in% (IN_CREP %>% pull(GW_wdid)) + +well_df <- well_df %>% group_by_at(c(c("GW_wdid","Year","AF"),SW_COL)) %>% summarize(IN_CREP=max(IN_CREP),CREP_Year=min(CREP_Year,na.rm=TRUE),CROP_AREA=sum(CROP_AREA),POT=sum(POT),ALF=sum(ALF),GP=sum(GP),SG=sum(SG),OT=sum(OT),SBD1=max(SBD1),SBD2=max(SBD2),SBD3=max(SBD3),SBD4=max(SBD4),SBD5=max(SBD5),SBD6=max(SBD6)) +well_df$CREP_Year[is.infinite(well_df$CREP_Year)] <- NA +temp <- well_df %>% group_by(GW_wdid) %>% mutate(ST_YEAR=min(Year)) +well_df <- well_df %>% select(GW_wdid,IN_CREP,CREP_Year,SBD1,Year,AF,CROP_AREA,POT,ALF,SG,GP,OT,all_of(SW_COL),SBD2,SBD3,SBD4,SBD5,SBD6) + +LIST2009 <- well_df %>% ungroup %>% filter(Year==2009) %>% pull(GW_wdid) %>% unique +MISSING2009 <- well_df[which(!(well_df$GW_wdid %in% LIST2009)),] %>% ungroup +MISSING2009 <- MISSING2009 %>% select(GW_wdid,IN_CREP,CREP_Year,SBD1,SBD2,SBD3,SBD4,SBD5,SBD6) %>% unique %>% mutate(Year=2009,AF=0,CROP_AREA=0,POT=0,ALF=0,SG=0,GP=0,OT=0) %>% left_join(WSW) +well_df <- rbind(MISSING2009,well_df) + +for(x in 2009:2021){ + S_YEAR <- x; + MISS <- well_df %>% group_by(GW_wdid,IN_CREP,CREP_Year,SBD1,SBD2,SBD3,SBD4,SBD5,SBD6) %>% filter(min(Year)<=S_YEAR,!any(Year==S_YEAR+1)) %>% select(GW_wdid,IN_CREP,CREP_Year,SBD1,SBD2,SBD3,SBD4,SBD5,SBD6) %>% ungroup %>% unique %>% mutate(Year=S_YEAR+1,AF=0,CROP_AREA=0,POT=0,ALF=0,SG=0,GP=0,OT=0) %>% left_join(WSW) + + well_df <- rbind(well_df,MISS) + print("done") +} +well_df <- well_df %>% group_by(GW_wdid,Year) %>% mutate(AF=sum(AF),IN_CREP=max(IN_CREP),CREP_Year=min(CREP_Year),CROP_AREA=sum(CROP_AREA),BOTH=ifelse(n()>1,1,0),POT=sum(POT),ALF=sum(ALF),SG=sum(SG),GP=sum(GP),OT=sum(OT),SBD1=max(SBD1),SW_2000812=max(SW_2000812),SW_2000623=max(SW_2000623),SW_2000631=max(SW_2000631),SW_2000753=max(SW_2000753),SW_2200593=max(SW_2200593),SW_3500570=max(SW_3500570),SW_2000829=max(SW_2000829),SW_2000798=max(SW_2000798),SW_2100521=max(SW_2100521),SW_2200541=max(SW_2200541),SBD1=max(SBD1),SBD2=max(SBD2),SBD3=max(SBD3),SBD4=max(SBD4),SBD5=max(SBD5),SBD6=max(SBD6)) %>% unique %>% ungroup + +well_df$Post_11 <- ifelse(well_df$Year>=2011,1,0) +well_df <- well_df %>% mutate(POST=ifelse(is.na(CREP_Year),0,ifelse(Year >=CREP_Year,1,0))) %>% select(GW_wdid,Year,POST,everything()) +well_df <- well_df %>% group_by(GW_wdid) %>% mutate(M_AF=max(AF),M_CROP=max(CROP_AREA)) %>% ungroup +well_df <- well_df %>% filter(M_AF!=0) + #only 2 2005 entries +well_df <- well_df %>% filter(Year>2005) + +##############Well Distances +AF_WELL <- well_df %>% group_by(GW_wdid) %>% summarize(AF=max(AF)) %>% arrange(GW_wdid) +well_list <- AF_WELL %>% filter(AF>0) %>% pull(GW_wdid) +WD <- read_csv("Data/Input_Data/Wells_Distance_Matrix.csv",col_types = cols(ID=col_character())) %>% rename(GW_wdid=ID) +##Keep only well records with pumping data, then arrange by well number for consistency +NM_SAVE <- WD[WD$GW_wdid %in% t(well_list),1] +###DISTANCE MATRIX +WD <- WD[WD$GW_wdid %in% well_list,colnames(WD) %in% well_list] +WD_FULL <- cbind(NM_SAVE,WD) %>% as_tibble + +##Data Clean +well_df$Yearf <- as.factor(well_df$Year) +well_df$GW_wdid <- as.factor(well_df$GW_wdid) +well_df$LAF <- log(well_df$AF+0.00001) +well_df <- well_df %>% rename(CREP_YEAR=CREP_Year) +well_df <- well_df %>% filter(Year!=2022) +well_df <- well_df %>% mutate(CREP_rel_year=Year-CREP_YEAR,CREP_rel_year=replace_na(CREP_rel_year,Inf)) +###Logit model of Fallowing +well_df$CREP_PRE <- ifelse(well_df$CREP_rel_year %in% c(-1,2),1,0) +well_df$CREP_PRE1 <- ifelse(well_df$CREP_rel_year %in% c(-1),1,0) +well_df$CREP_PRE2 <- ifelse(well_df$CREP_rel_year %in% c(-2),1,0) +well_df$CREP_PRE3 <- ifelse(well_df$CREP_rel_year %in% c(-3),1,0) +well_df$SUNAB_CREP <- ifelse(well_df$IN_CREP==0,100000,well_df$CREP_YEAR) +well_df <- well_df %>% left_join(IN_CREP %>% group_by(GW_wdid) %>% summarize(PERM_CREP=max(PERM_CREP),TEMP_CREP=as.numeric(min(PERM_CREP)==0))) %>% mutate(PERM_CREP=ifelse(is.na(PERM_CREP),0,PERM_CREP),TEMP_CREP=ifelse(is.na(TEMP_CREP),0,TEMP_CREP)) +PERM_CREP <- IN_CREP %>% filter(PERM_CREP==1) %>% group_by(GW_wdid) %>% summarize(PERM_YEAR=min(CREP_Year)) +TEMP_CREP <- IN_CREP %>% filter(PERM_CREP==0) %>% group_by(GW_wdid) %>% summarize(TEMP_YEAR=min(CREP_Year)) +well_df <- well_df %>% left_join(TEMP_CREP ) %>% mutate(ifelse(is.na(TEMP_CREP),0,TEMP_CREP)) +well_df <- well_df %>% left_join(PERM_CREP ) %>% mutate(ifelse(is.na(PERM_CREP),0,PERM_CREP)) +well_df <- well_df %>% mutate(CREP_YEAR=ifelse(is.na(CREP_YEAR),Inf,CREP_YEAR),TEMP_YEAR=ifelse(is.na(TEMP_YEAR),Inf,TEMP_YEAR),PERM_YEAR=ifelse(is.na(PERM_YEAR),Inf,PERM_YEAR)) +well_df$pumping_fee <- ifelse(well_df$Year==2011,45,NA) +well_df$pumping_fee <- ifelse(well_df$Year>2011 & well_df$Year<2018,75,well_df$pumping_fee) +well_df$pumping_fee <- ifelse(well_df$Year==2018 | well_df$Year==2019,90,well_df$pumping_fee) +well_df$pumping_fee <- ifelse(well_df$Year>2019 ,150,well_df$pumping_fee) +well_df$Post_16 <- ifelse(well_df$Year>=2016,1,0) +well_df$Post_11 <- ifelse(well_df$Year>=2011,1,0) + +#df_half$SUNAB <- ifelse(df_half$CLOSE_CREP==0,100000,df_half$CREP_CLOSE_START_YEAR) +well_df$FALLOW <- ifelse(well_df$AF==0,1,0) +well_df$BETWEEN<- ifelse(well_df$IN_CREP==1 & well_df$Year% left_join(FALL_PROG %>% mutate(GW_wdid=as.character(GW_wdid))) %>% rename(IN_FALL_PROG=FALL_PROG)%>% mutate(IN_FALL_PROG=ifelse(is.na(IN_FALL_PROG),0,1),FALL_PROG=ifelse(IN_FALL_PROG==1 & Year>=FALL_PROG_YEAR,1,0 )) +well_df$SBD1_Post <- ifelse(well_df$Post_11 & well_df$SBD1==1,1,0) +well_df$CREP_Post <- ifelse(well_df$Post_11 & well_df$IN_CREP==1,1,0) +well_df$SBD_ref_11 <- ifelse(well_df$SBD1==1,well_df$Year,Inf) +well_df$CREP_ref_11 <- ifelse(well_df$IN_CREP==1 & well_df$POST!=1,well_df$Year,Inf) +well_df$IND_CREP11 <- i(well_df$CREP_ref_11,ref=c(2009,2010,Inf)) +well_df$IND_SBD11 <- i(well_df$SBD_ref_11,ref=c(2009,2010,Inf)) +well_df <- well_df %>% select(-CREP_ref_11,-SBD_ref_11) +well_df$FALL_PROG_YEAR <- ifelse(is.na(well_df$FALL_PROG_YEAR),Inf,well_df$FALL_PROG_YEAR) + +source("Scripts/Distance_Functions.r") +GET_DIS_DATA_SET <- function(DIS,WELL_DF=well_df){ + DIS_DF <- WELL_DF %>% left_join(GET_CREP_DAT(DIS) ) + DIS_DF$CLOSE_CREP <- ifelse(DIS_DF$CREP_CLOSE_YEAR<100000,1,0) + DIS_DF$CLOSE_BETWEEN <- ifelse(DIS_DF$Post_11 & DIS_DF$CLOSE_CREP==1 & DIS_DF$CREP_CLOSE_YEAR=DIS_DF$Year ,1,0) + #All Fallow plots + DIS_DF <- DIS_DF %>% left_join(DF_CLOSE_FALL(DIS)) + #Data from 4 year Fallow program of SBD1 + DIS_DF <- DIS_DF %>% left_join(GET_FALLOW_PROG_DAT(DIS)) + DIS_DF <- DIS_DF %>% left_join(GET_FALLOW_PROG_DAT(DIS)) + DIS_DF$FALL_PROG_CLOSE_YEAR <- ifelse(is.na(DIS_DF$FALL_PROG_CLOSE_YEAR),Inf,DIS_DF$FALL_PROG_CLOSE_YEAR) + DIS_DF$FALL_PROG_CLOSE_POST <- ifelse(DIS_DF$FALL_PROG==1 & DIS_DF$FALL_PROG_CLOSE_YEAR>=DIS_DF$Year,1,0) + return(DIS_DF) +} +df_quart <- GET_DIS_DATA_SET(5280/4) +df_half <- GET_DIS_DATA_SET(5280/2) +df_mile <- GET_DIS_DATA_SET(5280) +df_two <- GET_DIS_DATA_SET(2*5280) +#########################ADDING IN PROBIT MODEL OF ENTERING CREP DUE TO FEE +PRE_DF <- df_half %>% group_by(GW_wdid,IN_CREP,CLOSE_CREP)%>% filter(Year<2011,SBD1==1) %>% summarize(M_AF=max(AF),PRE_AF=mean(AF),PRE_CROP=mean(POT+ALF+SG+OT),PRE_POT=mean(POT)/PRE_CROP,PRE_ALF=mean(ALF)/PRE_CROP,PRE_SG=mean(SG)/PRE_CROP,PRE_OT=mean(OT)/PRE_CROP) %>% filter(PRE_CROP>0) +POST_DF <- df_half %>% group_by(GW_wdid,IN_CREP,CLOSE_CREP)%>% filter(Year>=2011,Year<2014,SBD1==1) %>% summarize(POST_AF=mean(AF),POST_CROP=mean(POT+ALF+SG+OT,na.rm=TRUE),POST_POT=mean(POT,na.rm=TRUE)/POST_CROP,POST_ALF=mean(ALF,na.rm=TRUE)/POST_CROP,POST_SG=mean(SG,na.rm=TRUE)/POST_CROP,POST_OT=mean(OT,na.rm=TRUE)/POST_CROP) +SW_DF <- df_half[,c(1,19:28)] %>% unique +logit_df <- PRE_DF %>% left_join(POST_DF) %>% mutate(DIFF_AF=POST_AF-PRE_AF,DIFF_CROP=PRE_CROP-POST_CROP) %>% left_join(SW_DF) diff --git a/Hydrobase.r b/2_Hydrobase.r old mode 100644 new mode 100755 similarity index 100% rename from Hydrobase.r rename to 2_Hydrobase.r diff --git a/Data/Readme_Data.md b/Data/Readme_Data.md new file mode 100644 index 0000000..5fc5aab --- /dev/null +++ b/Data/Readme_Data.md @@ -0,0 +1,5 @@ +If no data is availble download the input files from my personal cloud drive. + + https://u.pcloud.link/publink/show?code=XZfpNd5ZUbz9DW46DYY6g31ADh41B5U089yy + +Download and extract this folder in the "Data" directory. diff --git a/Scripts/Distance_Functions.r b/Scripts/Distance_Functions.r new file mode 100644 index 0000000..562bf02 --- /dev/null +++ b/Scripts/Distance_Functions.r @@ -0,0 +1,75 @@ +WD <- as.matrix(WD) +GET_CUTOFF <- function(DIST){ + WD_MAT <- as.matrix(WD) + WD_PICK <- as.matrix(WD) + WD_PICK[WD_MAT<=DIST] <- TRUE + WD_PICK[WD_MAT>DIST] <- FALSE + diag(WD_PICK) <- FALSE + return(WD_PICK) +} +COND_MAT <- function(SEL_LIST,DIST){ + NUM_RES <- colSums(GET_CUTOFF(DIST)*(colnames(WD) %in% SEL_LIST)) + return(cbind(colnames(WD),NUM_RES) %>% as_tibble%>% mutate(NUM_RES=as.numeric(NUM_RES),IN=ifelse(NUM_RES>0,TRUE,FALSE)) %>% rename(GW_wdid=V1)) +} + + +GET_CREP_DAT <- function(DIS){ + for(i in c(2014,2015,2016,2017,2018,2019,2020)){ + SEL_LIST <- IN_CREP %>% filter(CREP_Year==i,PERM_CREP==0) %>% pull(GW_wdid) + if(!exists("RES")){RES <- COND_MAT(SEL_LIST,DIS)%>% mutate(CREP_CLOSE_YEAR=i)}else{RES <- rbind(RES,COND_MAT(SEL_LIST,DIS)%>% mutate(CREP_CLOSE_YEAR=i))} + } + RES$PERM_CREP <- 0 + for(i in c(2014,2015,2016,2018,2020,2021)){ + SEL_LIST <- IN_CREP %>% filter(CREP_Year==i,PERM_CREP==1) %>% pull(GW_wdid) + RES <- rbind(RES,COND_MAT(SEL_LIST,DIS)%>% mutate(CREP_CLOSE_YEAR=i,PERM_CREP=1)) + } + for(i in 2014:2021){ + CREP <- RES %>% filter(CREP_CLOSE_YEAR<=i) %>% group_by(GW_wdid) %>% summarize(CREP_CLOSE_NUM=sum(NUM_RES),CREP_CLOSE_TREAT=max(IN),Year=i) + TEMP <- RES %>% filter(CREP_CLOSE_YEAR<=i,PERM_CREP==0) %>% group_by(GW_wdid) %>% summarize(TEMP_CLOSE_NUM=sum(NUM_RES),TEMP_CLOSE_TREAT=max(IN),Year=i) + PERM <- RES %>% filter(CREP_CLOSE_YEAR<=i,PERM_CREP==1) %>% group_by(GW_wdid) %>% summarize(PERM_CLOSE_NUM=sum(NUM_RES),PERM_CLOSE_TREAT=max(IN),Year=i) + + ALL <- full_join(full_join(CREP,TEMP),PERM) + if(!exists("NUM_RES")){NUM_RES <- ALL}else{NUM_RES <- rbind(NUM_RES,ALL)} + } + for(i in 2009:2013){ + temp <- NUM_RES %>% filter(Year==2014) + temp[,-1] <- 0 + temp$Year <- i + NUM_RES <- rbind(NUM_RES,temp) + } + + POLICY_SUM <- full_join(RES %>% filter(PERM_CREP==0) %>% group_by(GW_wdid) %>% summarize(TEMP_CLOSE_YEAR=min(ifelse(IN,CREP_CLOSE_YEAR,Inf))),RES %>% filter(PERM_CREP==1) %>% group_by(GW_wdid) %>% summarize(PERM_CLOSE_YEAR=min(ifelse(IN,CREP_CLOSE_YEAR,Inf)))) + POLICY_SUM <- POLICY_SUM %>% left_join(RES %>% group_by(GW_wdid) %>% summarize(CREP_CLOSE_YEAR=min(ifelse(IN,CREP_CLOSE_YEAR,Inf)))) + + POLICY_SUM <- POLICY_SUM %>% full_join(NUM_RES) %>% select(GW_wdid,Year,everything()) %>% arrange(GW_wdid,Year) + return(POLICY_SUM) +} +##### +GET_CLOSE_FALLOW <- function(YEAR,DIST){ + OR <- colnames(WD) + FALL <- well_df %>% filter(Year==YEAR) %>% arrange(GW_wdid==OR) %>% mutate(FALLOW=(AF==0)) %>% pull(FALLOW) + WD <- as.matrix(WD) + CLOSE_FALLOW <- cbind(OR,rep(YEAR,nrow(WD)),rowSums(WD[,FALL] % as_tibble + colnames(CLOSE_FALLOW) <- c("GW_wdid","Year","NUM_CLOSE_FALLOW") + CLOSE_FALLOW$Year <- as.numeric(CLOSE_FALLOW$Year) + CLOSE_FALLOW$NUM_CLOSE_FALLOW<- as.integer(CLOSE_FALLOW$NUM_CLOSE_FALLOW) + return(CLOSE_FALLOW) +} +DF_CLOSE_FALL <- function(DIST){ + CLOSE_FAL <- GET_CLOSE_FALLOW(2009,DIST) + for(x in 2010:2021){ + CLOSE_FAL <- rbind(CLOSE_FAL,GET_CLOSE_FALLOW(x,DIST)) + } + CLOSE_FAL <- CLOSE_FAL %>% mutate(HAS_CLOSE_FALL=ifelse(NUM_CLOSE_FALLOW>0,1,0)) + return(CLOSE_FAL) +} +################################################################################### +GET_FALLOW_PROG_DAT <- function(DIS){ + for(i in c(2020,2021)){ + SEL_LIST <- FALL_PROG %>% filter(FALL_PROG_YEAR==i) %>% pull(GW_wdid) + if(!exists("RES")){RES <- COND_MAT(SEL_LIST,DIS)%>% mutate(FALL_PROG_CLOSE_YEAR=i)}else{RES <- rbind(RES,COND_MAT(SEL_LIST,DIS)%>% mutate(FALL_PROG_CLOSE_YEAR=i))} + } +RES <- RES %>% filter(IN) %>% select(GW_wdid,FALL_PROG_CLOSE_YEAR) +RES <- RES %>% group_by(GW_wdid) %>% mutate(FALL_PROG_CLOSE_YEAR=min(FALL_PROG_CLOSE_YEAR)) %>% unique +return(RES) +}