Update data to better observe SBD1 entrance time

1_Create_CREP_Data.r

@@ -135,3 +135,6 @@ ALL_PROGRAMS <- ALL_PROGRAMS%>% pivot_wider(values_from=active,names_from=progra
 ALL_PROGRAMS <- ALL_PROGRAMS %>% mutate(CREP_any=ifelse(CREP_perm+CREP_temp>0,1,0)) %>% select(wdid,year,CREP_any,everything())
 
 write_csv(ALL_PROGRAMS,"Data/Output_Data/Fallow_Program_Data.csv")
+saveRDS(ALL_PROGRAMS,"Data/Output_Data/Fallow_Program_Data.rds")
+
+

1_Crop_Choice.r

@@ -1,34 +0,0 @@
-library(tidyverse)
-CROPS <- rbind(read_csv("Data/Crop_Choice/IRRIG_2002.csv"),read_csv("Data/Crop_Choice/IRRIG_2005.csv"))
-DITCH_PARCEL <- CROPS[,c(4,11:19)] %>% pivot_longer(-PARCEL_ID,values_to="ditch") %>% filter(!is.na(ditch)) %>% select(-name)
-WELL_PARCEL <- 	CROPS[,c(4,seq(21,80,by=3))] %>% pivot_longer(-PARCEL_ID,values_to="wdid") %>% filter(!is.na(wdid)) %>% select(-name)
-WELL_DITCHES <- WELL_PARCEL %>% left_join(DITCH_PARCEL)  %>% mutate(ditch=as.character(ditch))
-NAMED_DITCHES <- (WELL_DITCHES  %>% filter(!is.na(ditch))   %>% group_by(ditch)  %>% summarize(size=n()) %>% arrange(desc(size)) )[1:10,] %>% pull(ditch)
-WELL_DITCHES[!(WELL_DITCHES$ditch  %in% NAMED_DITCHES) & !is.na(WELL_DITCHES$ditch) ,3] <- "Other"
-WELL_DITCHES <- WELL_DITCHES[-which(duplicated(WELL_DITCHES)),]
-WELL_DITCHES <- WELL_DITCHES  %>% filter(!is.na(ditch),!is.na(wdid)) 
-WELL_DITCHES$ROW_NUM <- 1:nrow(WELL_DITCHES)
-WELL_DITCHES <- WELL_DITCHES %>% select(-PARCEL_ID) %>% unique %>% pivot_wider(values_from=ditch,names_from=ditch,names_prefix="ditch_")
-WELL_DITCHES[,-1:-2][!is.na(WELL_DITCHES[,-1:-2])] <- '1'
-WELL_DITCHES[,-1:-2][is.na(WELL_DITCHES[,-1:-2])] <- '0'
-WELL_DITCHES <-  WELL_DITCHES %>% mutate_if(is.character,as.numeric) %>% select(-ROW_NUM)  %>% unique
-WELL_DITCHES <- WELL_DITCHES   %>% group_by(wdid) %>%   mutate(across(colnames(WELL_DITCHES)[2:12], max, na.rm = TRUE)) %>% ungroup  %>% unique
-dir.create("Data/Output_Data",showWarnings=FALSE)
-write_csv(WELL_DITCHES,"Data/Output_Data/Ditch_Indicators.csv")
-###########Determine what percentage of land supplied by a well goes to each crop, across 2002 and 2005 
-CROP <- CROPS %>% select(CAL_YEAR,PARCEL_ID,ACRES,CROP_TYPE)
-CROP$CROP_TYPE <- ifelse(CROP$CROP_TYPE=='WHEAT_FALL','SMALL_GRAINS',CROP$CROP_TYPE)
-CROP$CROP_TYPE <- ifelse(CROP$CROP_TYPE=='NEW_ALFALFA','ALFALFA',CROP$CROP_TYPE)
-CROP$CROP_TYPE <- ifelse(CROP$CROP_TYPE %in% c('COVER_CROP','VEGETABLES'),'OTHER',CROP$CROP_TYPE)
-
-CROP_2002 <- WELL_PARCEL %>% left_join(CROP %>% filter(CAL_YEAR==2002))%>% filter(!is.na(CROP_TYPE)) %>% unique 
-CROP_2005 <- WELL_PARCEL %>% left_join(CROP %>% filter(CAL_YEAR==2005))%>% filter(!is.na(CROP_TYPE)) %>% unique 
-
-CROP_2002 <- CROP_2002 %>% group_by(wdid,CROP_TYPE,CAL_YEAR) %>% summarize(ACRES=sum(ACRES)) %>% group_by(wdid,CAL_YEAR) %>% mutate(PERCENT=ACRES/sum(ACRES))  %>% arrange(wdid) %>% ungroup %>% select(-ACRES) %>% pivot_wider(values_from=PERCENT,names_from=CROP_TYPE,names_prefix="PER_") %>% replace(is.na(.), 0)
-CROP_2005 <- CROP_2005 %>% group_by(wdid,CROP_TYPE,CAL_YEAR) %>% summarize(ACRES=sum(ACRES)) %>% group_by(wdid,CAL_YEAR) %>% mutate(PERCENT=ACRES/sum(ACRES))  %>% arrange(wdid) %>% ungroup %>% select(-ACRES) %>% pivot_wider(values_from=PERCENT,names_from=CROP_TYPE,names_prefix="PER_") %>% replace(is.na(.), 0)
-PRE_2009_CROPS <- full_join(CROP_2002 %>% select(-CAL_YEAR) ,CROP_2005 %>% select(-CAL_YEAR)) %>% group_by(wdid) %>% summarize(across(colnames(CROP_2002)[3:7], mean, na.rm = TRUE))
-PRE_2009_CROPS <- PRE_2009_CROPS  %>% clean_names()
-write_csv(PRE_2009_CROPS ,"Data/Output_Data/Crops_Before_2009.csv")
-
-
-

2_Proc.r

@@ -1,11 +1,7 @@
-#install.packages("devtools")
-#devtools::install_github("anguswg-ucsb/cdssr")
-#library(cdssr)
 library(tidyverse)
 library(janitor)
-#install.packages("janitor")
+
-#API <- 'PwevUQJCStcYZfqrOYbuyztmNPlUJWby'
+#Read the data files which include the current list of wells by subdistrict, to create indicator variables. This is the list of wells currently in the subdistrict using hydrobase search filter. Note that the time of entry is not found here and added later.
-	read_csv("Data/SBD_Data/StructureList_SBD2.csv")
 SBD_LINK <- rbind(read_csv("Data/SBD_Data/StructureList_SBD1.csv") %>% mutate(SBD=1),
 read_csv("Data/SBD_Data/StructureList_SBD2.csv") %>% mutate(SBD=2),
 read_csv("Data/SBD_Data/StructureList_SBD3.csv") %>% mutate(SBD=3),
@@ -14,21 +10,33 @@ read_csv("Data/SBD_Data/StructureList_SBD5.csv") %>% mutate(SBD=5),
 read_csv("Data/SBD_Data/StructureList_SBD6.csv") %>% mutate(SBD=6)) %>% clean_names %>% mutate(wdid=as.character(wdid))
 SBD_LINK <- SBD_LINK %>% pivot_wider(values_from=sbd,names_from=sbd,names_prefix="SBD") %>% replace(is.na(.), 0)
 SBD_LINK[2:7][SBD_LINK[2:7]>0] <- 1
+SBD1_WELLS <- SBD_LINK %>% select(wdid)  %>% unique
+############Data from each of the Annuall Replacement Plans or in the case of 2000 the Plan of Water Management. A list of wells in Subdistrict one was collected to create a time depdent indicator for being in SBD1.
+SBD1_WELLS <- read_csv("Data/ARP/SBD1_Wells.csv")%>% group_by(wdid) %>% filter(year==min(year)) %>% mutate(wdid=as.character(wdid))  %>% rename(SBD1_year=year) %>% ungroup %>% inner_join(SBD1_WELLS)
+SBD1_WELLS %>% group_by(SBD1_year) %>% summarize(n())
+SBD1_WELLS %>% filter(wdid=='2706148')
+######################Static Data
+WELL_DATA <- read_csv("Data/Structure_Data/Structures_with_Diversions.csv")%>% clean_names() %>% mutate(wdid=as.character(wdid)) %>% select(wdid,contacts,latitude,longitude) #Start with well data
 
-#################
+STATIC_DATA <- WELL_DATA%>% left_join(SBD1_WELLS) %>% left_join(SBD_LINK) %>% replace(is.na(.), 0) #Include subdistrict indicators
-#STRUCTURES <- 'https://dwr.state.co.us/Rest/GET/api/v2/structures/?format=csv&fields=wdid%2CciuCode%2CstructureType&division=3&pageSize=500000&apiKey=PwevUQJCStcYZfqrOYbuyztmNPlUJWby'
+STATIC_DATA$SBD1_year <- ifelse(STATIC_DATA$SBD1_year==0,Inf,STATIC_DATA$SBD1_year  )
-#STRUCTURES<- read_csv(STRUCTURES,skip=2 )
-#WELLS <- STRUCTURES %>% filter(structuretype=='WELL')
-
-WELL_DATA <- read_csv("Data/Structure_Data/Structures_with_Diversions.csv")%>% clean_names()  %>% mutate(wdid=as.character(wdid)) %>% select(wdid,contacts) #Start with well data
-STATIC_DATA <- WELL_DATA %>% left_join(SBD_LINK) %>% replace(is.na(.), 0) #Include subdistrict indicators
 STATIC_DATA  <- STATIC_DATA  %>% left_join(read_csv("Data/Output_Data/Ditch_Indicators.csv")%>% mutate(wdid=as.character(wdid)))%>% replace(is.na(.), 0)
 
 STATIC_DATA <- STATIC_DATA%>% left_join(read_csv("Data/Output_Data/Crops_Before_2009.csv") %>% mutate(wdid=as.character(wdid))) #Add crop data
 STATIC_DATA$CROPS_PRE_2009 <- ifelse(is.na(STATIC_DATA$per_alfalfa),0,1) #Make an indicator to tell if crops were grown in 2002 or 2005, or if not crop data was available.
 STATIC_DATA <- STATIC_DATA %>% replace(is.na(.), 0)
-
+SBD1_ORIG_WELLS <- STATIC_DATA %>% filter(SBD1_year==2009) %>% select(wdid,latitude,longitude)
 write_csv(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.csv")
+write_csv(SBD1_ORIG_WELLS,"Data/Output_Data/Original_SBD1_Well_Location.csv")
+#######Original Wells based on map
+	#This data comes from a manual map of SBD1 in 2009 made in QGIS. Wells are intersected with the hand drawn layer. The wells identified as being in the SBD1 region at this time should be assumed to have entered the subdistrict at the start of the program (2006 with first pumping data in 2009). This fixes the problem that some wells are identified as being in SBD1 in ARP reports but only because they were missed in earlier reports. It looks like most updated wells were missed due to category issues, such as not knowing the well was used for agriculture. This process should find most wells that were enrolled at the start, and later picked up so as to treat them as 2006/2009 wells instead of treating them as having entered in say 2016 or 2020 when they were first picked up by the SBD1 list. 
+SBD1_FROM_MAP <- read_csv("Data/SBD_Data/Intersect_SBD1_Area_Wells.csv")
+SBD1_FROM_MAP <-  SBD1_FROM_MAP %>% pull(WELL_ID) %>% unique 
+SBD1_FROM_MAP <- SBD1_FROM_MAP[(SBD1_FROM_MAP %in% (STATIC_DATA %>% filter(SBD1==1))$wdid  )] #Remove any wells that are still not in SBD1. These were most likely picked up in error. For example if part of the parcel is in SBD1 but the well serves another parcel outside of the area, or because I drew the boundary too large.
+
+STATIC_DATA[STATIC_DATA$wdid %in% SBD1_FROM_MAP ,"SBD1_year"] <- 2009
+write_csv(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.csv")
+saveRDS(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.rds")
 
 
 ###########Collect Pumping Data
@@ -37,8 +45,6 @@ ALL <- length(LIST )
 
 if(file.exists("./Data/Output_Data/Div3_Pumping_Data.csv")){LIST <-LIST[-which(LIST %in%	t(read.csv("./Data/Output_Data/Div3_Pumping_Data.csv")[,1] %>% unique))]}
 CURRENT <- length(LIST )
-ALL
-CURRENT 
 
 WITH_API_KEY <- FALSE 
 for(C_WELL in LIST ){
@@ -52,12 +58,9 @@ PUMPING <- read_csv("./Data/Output_Data/Div3_Pumping_Data.csv") %>% unique %>% r
 
 if(colnames(PUMPING)[1]!='wdid'){PUMPING <- read_csv("./Data/Output_Data/Div3_Pumping_Data.csv",col_names=c("wdid","year","AF")) %>% unique} #See if the names were already added, or if they need renamed
 POST_LIST <- PUMPING %>% pull(wdid) %>% unique
-PUMPING %>% group_by(wdid,year) %>% filter(n()>1) %>% arrange(wdid,year)
-POST_LIST[1]
 PUMPING$ROW <- 1:nrow(PUMPING)
 PUMPING <- PUMPING %>% pivot_wider(values_from=AF,names_from=year)%>% group_by(wdid) %>% summarize(across(as.character(2009:2025),\(x) mean(x, na.rm = TRUE))) %>% pivot_longer(-wdid,names_to='year',values_to='AF') %>% mutate(wdid=as.character(wdid)) %>% unique #Pivot to add zeros when a year is missing data
 write_csv(PUMPING,file="./Data/Output_Data/Div3_Pumping_Data.csv")
 ALL_DATA <- PUMPING %>% left_join(STATIC_DATA) %>% clean_names()
 write_csv(ALL_DATA,file="./Data/Output_Data/Full_Data_Set.csv")
 
-

3_Proc_Dynamic_Data.r

@@ -0,0 +1,41 @@
+library(tidyverse)
+library(janitor)
+
+###########Collect Pumping Data
+WELL_DATA <- read_csv("Data/Structure_Data/Structures_with_Diversions.csv")%>% clean_names() %>% mutate(wdid=as.character(wdid)) %>% select(wdid,contacts,latitude,longitude) #Start with well data
+STATIC_DATA <- readRDS("Data/Output_Data/Well_Level_Static_Data.rds")
+FALLOW_PROGRAM_DATA <- readRDS("Data/Output_Data/Fallow_Program_Data.rds")
+
+
+LIST <- WELL_DATA %>% pull(wdid) %>% unique 
+ALL <- length(LIST )
+
+if(file.exists("./Data/Output_Data/Div3_Pumping_Data.csv")){LIST <-LIST[-which(LIST %in%	t(read.csv("./Data/Output_Data/Div3_Pumping_Data.csv")[,1] %>% unique))]}
+CURRENT <- length(LIST )
+
+WITH_API_KEY <- FALSE 
+for(C_WELL in LIST ){
+	if(WITH_API_KEY){try(C_DAT <-	read_csv(paste0('https://dwr.state.co.us/Rest/GET/api/v2/structures/divrec/divrecyear/?format=csv&dateFormat=dateOnly&fields=wdid%2CdataMeasDate%2CdataValue&measUnits=ACFT&wcIdentifier=*Total+(Diversion)*&wdid=',C_WELL,'&apiKey=PwevUQJCStcYZfqrOYbuyztmNPlUJWby'),skip=2) %>% select(wdid,year=datameasuate,pumping=dataValue))} else{
+try(C_DAT <-	read_csv(paste0('https://dwr.state.co.us/Rest/GET/api/v2/structures/divrec/divrecyear/?format=csv&dateFormat=dateOnly&fields=wdid%2CdataMeasDate%2CdataValue&measUnits=ACFT&wcIdentifier=*Total+(Diversion)*&wdid=',C_WELL),skip=2) %>% select(wdid,year=dataMeasDate,pumping=dataValue))
+	}
+if(exists("C_DAT")){	if(C_DAT[[1,1]]==C_WELL){write_csv(C_DAT,append=TRUE,file="./Data/Output_Data/Div3_Pumping_Data.csv")}}
+	#closeAllConnections()
+}
+PUMPING <- read_csv("./Data/Output_Data/Div3_Pumping_Data.csv") %>% unique %>% replace(is.na(.), 0)
+
+if(colnames(PUMPING)[1]!='wdid'){PUMPING <- read_csv("./Data/Output_Data/Div3_Pumping_Data.csv",col_names=c("wdid","year","AF")) %>% unique} #See if the names were already added, or if they need renamed
+POST_LIST <- PUMPING %>% pull(wdid) %>% unique
+PUMPING$ROW <- 1:nrow(PUMPING)
+PUMPING <- PUMPING %>% pivot_wider(values_from=AF,names_from=year)%>% group_by(wdid) %>% summarize(across(as.character(2009:2025),\(x) mean(x, na.rm = TRUE))) %>% pivot_longer(-wdid,names_to='year',values_to='AF') %>% mutate(wdid=as.character(wdid)) %>% unique #Pivot to add zeros when a year is missing data
+#One value is slightly negative
+PUMPING$AF <- ifelse(PUMPING$AF<0,0,PUMPING$AF)
+write_csv(PUMPING,file="./Data/Output_Data/Div3_Pumping_Data.csv")
+
+ALL_DATA <- PUMPING %>% left_join(STATIC_DATA) %>% clean_names()
+DYNAMIC_DATA <- PUMPING %>% left_join(FALLOW_PROGRAM_DATA)%>% replace(is.na(.), 0)
+ALL_DATA <- DYNAMIC_DATA %>% left_join(STATIC_DATA) 
+
+write_csv(ALL_DATA,file="./Data/Output_Data/Full_Data_Set.csv")
+saveRDS(ALL_DATA,file="./Data/Output_Data/Full_Data_Set.rds")
+
+

4_Analysis.r

@@ -1,14 +1,29 @@
 library(tidyverse)
 library(fixest)
-WELLS <- read_csv("Data/Output_Data/Full_Data_Set.csv")
+DF <- readRDS("Data/Output_Data/Full_Data_Set.rds")
-WELLS$post <- WELLS$year>=2011
+DF$POST <- ifelse(DF$year>=2011,1,0)
-WELLS$I <- WELLS$year-2011
+
-feols(af~sbd1*post|ditch_2000812+ ditch_2000631+ditch_other+ditch_2000829+ditch_2000798+ditch_2000816+ditch_2000753+ditch_2000623+ditch_2200627+contacts+year+sbd2+sbd3+sbd4+sbd5+sbd6,data=WELLS)
+REG_LN1 <- feols(log(AF+0.001)~POST*SBD1,cluster="wdid",data=DF)
-iplot(feols(log(af+0.0001)~i(I,sbd1)|ditch_2000812+ ditch_2000631+ditch_other+ditch_2000829+ditch_2000798+ditch_2000816+ditch_2000753+ditch_2000623+ditch_2200627+contacts+year+sbd2+sbd3+sbd4+sbd5+sbd6,data=WELLS))
+REG_AF1 <- feols(AF~POST*SBD1,cluster="wdid",data=DF)
-iplot(feols(af~i(year,sbd1,2010)+per_alfalfa+per_potatoes+per_alfalfa|year+sbd2+sbd3+sbd4+sbd5+sbd6+ditch_2000812^year+ditch_2000631^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year,data=WELLS))
+REG_AF1 <- feols(AF~POST*SBD1+CREP_any+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_potatoes+per_alfalfa|SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,cluster="wdid",data=DF)
+REG_LN1 <- feols(log(AF+0.001)~POST*SBD1+CREP_any+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_potatoes+per_alfalfa|SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,cluster="wdid",data=DF)
+
+DF$SBD1 <- ifelse(DF$year>=DF$SBD1_year,1,0)
+DF
+DF$CROPS_PRE_2009
+
+REG_AF2 <- feols(AF~POST*SBD1+CREP_any+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_potatoes+per_alfalfa|SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,cluster="wdid",data=DF)
+REG_LN2 <- feols(log(AF+0.001)~POST*SBD1+CREP_any+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_potatoes+per_alfalfa|SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,cluster="wdid",data=DF)
 ?i
+DF$I_VAL <- DF$SBD1_year-as.numeric(DF$year)
+summary(DF$I_VAL)
+
+iplot( feols(AF~i(I_VAL,ref=-1)+CREP_any+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_potatoes+per_alfalfa|SBD2+SBD3+SBD4+SBD5+SBD6+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,cluster="wdid",data=DF))
+DF$SBD1_year <- ifelse(DF$SBD1_year>10000,-Inf,DF$SBD1_year)
+DF <- DF %>% group_by(wdid) %>% mutate(SBD1=max(SBD1))
 
 
-WELLS <-WELLS %>% group_by(year,sbd1) %>% summarize(af=sum(af)) %>% ungroup
+
-ggplot(WELLS,aes(x=year,y=af,group=sbd1,color=(sbd1)))+geom_point()+geom_line()
+etable(REG_AF1,REG_AF2,headers=c("Static SBD1","Dynamic SBD1"),fixef.group=list("Other Subdistricts"=c("SBD"),"Ditch-Year"="ditch"),file="Reg.tex",export="Reg.png",replace=TRUE)
+coefplot(REG_AF1,REG_AF2)