From 6024203d1c5f381535be6b8f03eb3ce2def23399 Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Thu, 7 May 2026 12:55:23 -0600 Subject: [PATCH] Update data to better observe SBD1 entrance time --- 2_Create_CREP_Data.r => 1_Create_CREP_Data.r | 3 ++ 1_Crop_Choice.r | 34 ---------------- 3_Hydrobase.r => 2_Proc.r | 43 +++++++++++--------- 3_Proc_Dynamic_Data.r | 41 +++++++++++++++++++ 4_Analysis.r | 31 ++++++++++---- 5 files changed, 90 insertions(+), 62 deletions(-) rename 2_Create_CREP_Data.r => 1_Create_CREP_Data.r (99%) delete mode 100644 1_Crop_Choice.r rename 3_Hydrobase.r => 2_Proc.r (51%) create mode 100644 3_Proc_Dynamic_Data.r diff --git a/2_Create_CREP_Data.r b/1_Create_CREP_Data.r similarity index 99% rename from 2_Create_CREP_Data.r rename to 1_Create_CREP_Data.r index 36ae370..de616dd 100644 --- a/2_Create_CREP_Data.r +++ b/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") + + diff --git a/1_Crop_Choice.r b/1_Crop_Choice.r deleted file mode 100644 index efbd83d..0000000 --- a/1_Crop_Choice.r +++ /dev/null @@ -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") - - - diff --git a/3_Hydrobase.r b/2_Proc.r similarity index 51% rename from 3_Hydrobase.r rename to 2_Proc.r index a0d21d3..77f4aa2 100644 --- a/3_Hydrobase.r +++ b/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_csv("Data/SBD_Data/StructureList_SBD2.csv") + +#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. 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 -################# -#STRUCTURES <- 'https://dwr.state.co.us/Rest/GET/api/v2/structures/?format=csv&fields=wdid%2CciuCode%2CstructureType&division=3&pageSize=500000&apiKey=PwevUQJCStcYZfqrOYbuyztmNPlUJWby' -#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 <- WELL_DATA%>% left_join(SBD1_WELLS) %>% left_join(SBD_LINK) %>% replace(is.na(.), 0) #Include subdistrict indicators +STATIC_DATA$SBD1_year <- ifelse(STATIC_DATA$SBD1_year==0,Inf,STATIC_DATA$SBD1_year ) +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") - diff --git a/3_Proc_Dynamic_Data.r b/3_Proc_Dynamic_Data.r new file mode 100644 index 0000000..39db45c --- /dev/null +++ b/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") + + diff --git a/4_Analysis.r b/4_Analysis.r index ab6fde7..92cb000 100644 --- a/4_Analysis.r +++ b/4_Analysis.r @@ -1,14 +1,29 @@ library(tidyverse) library(fixest) -WELLS <- read_csv("Data/Output_Data/Full_Data_Set.csv") -WELLS$post <- WELLS$year>=2011 -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) -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)) -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)) +DF <- readRDS("Data/Output_Data/Full_Data_Set.rds") +DF$POST <- ifelse(DF$year>=2011,1,0) + +REG_LN1 <- feols(log(AF+0.001)~POST*SBD1,cluster="wdid",data=DF) +REG_AF1 <- feols(AF~POST*SBD1,cluster="wdid",data=DF) +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)