alex/CREP_Policy_Interaction_Analysis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added code to seperate water class for irrigation

Browse Source
This commit is contained in:
Alex Gebben Work 2026-05-14 16:56:39 -06:00
parent 52d3ea6ff0
commit f8fb0df3f3
2 changed files with 57 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
2_Proc.r
View File

@ -15,6 +15,14 @@ SBD1_WELLS <- SBD_LINK %>% select(wdid) %>% unique
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 <- 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 %>% group_by(SBD1_year) %>% summarize(n())
SBD1_WELLS %>% filter(wdid=='2706148') SBD1_WELLS %>% filter(wdid=='2706148')
#######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% SBD1_WELLS$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.
SBD1_WELLS[SBD1_WELLS$wdid %in% SBD1_FROM_MAP,"SBD1_year"] <- 2011
SBD1_WELLS %>% filter(SBD1_year!=2011)
######################Static Data ######################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 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
@ -23,22 +31,25 @@ STATIC_DATA$SBD1_year <- ifelse(STATIC_DATA$SBD1_year==0,Inf,STATIC_DATA$SBD1_ye
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/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 <- 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$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 no crop data was available.
STATIC_DATA <- STATIC_DATA %>% replace(is.na(.), 0) STATIC_DATA <- STATIC_DATA %>% replace(is.na(.), 0)
SBD1_ORIG_WELLS <- STATIC_DATA %>% filter(SBD1_year==2009) %>% select(wdid,latitude,longitude) ###Make sure all CREP wells are included in SBD1
write_csv(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.csv") PROGRAM_WELLS <- c(read_csv("Data/CREP_and_Fallow/Cleaned/Fallow_Program_Data/Bill_SB22_Fallow_Payments.csv")$wdid,
write_csv(SBD1_ORIG_WELLS,"Data/Output_Data/Original_SBD1_Well_Location.csv") read_csv("Data/CREP_and_Fallow/Cleaned/Fallow_Program_Data/CREP.csv")$wdid,
#######Original Wells based on map read_csv("Data/CREP_and_Fallow/Cleaned/Fallow_Program_Data/SBD1_Half_Fallow_Program.csv")$wdid,
#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. read_csv("Data/CREP_and_Fallow/Cleaned/Fallow_Program_Data/SBD1_Temporary_Fallow_Payments.csv")$wdid)
SBD1_FROM_MAP <- read_csv("Data/SBD_Data/Intersect_SBD1_Area_Wells.csv") STATIC_DATA[STATIC_DATA$wdid %in% PROGRAM_WELLS,"SBD1_year"] <- 2011
SBD1_FROM_MAP <- SBD1_FROM_MAP %>% pull(WELL_ID) %>% unique STATIC_DATA[STATIC_DATA$wdid %in% PROGRAM_WELLS,"SBD1"] <- 1
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. AUG <- read_csv("Data/Augmentation/Augmentation.csv") %>% mutate(wdid=as.character(wdid)) %>% select(wdid)
STATIC_DATA$augmentation_well <- 0
STATIC_DATA[STATIC_DATA$wdid %in% AUG, "augmentation_well" ] <- 1
STATIC_DATA[STATIC_DATA$wdid %in% SBD1_FROM_MAP ,"SBD1_year"] <- 2009
write_csv(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.csv") write_csv(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.csv")
saveRDS(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.rds") saveRDS(STATIC_DATA,"Data/Output_Data/Well_Level_Static_Data.rds")
###########Collect Pumping Data ###########Collect Pumping Data
LIST <- WELL_DATA %>% pull(wdid) %>% unique LIST <- WELL_DATA %>% pull(wdid) %>% unique
ALL <- length(LIST ) ALL <- length(LIST )
@ -48,8 +59,8 @@ CURRENT <- length(LIST )
WITH_API_KEY <- FALSE WITH_API_KEY <- FALSE
for(C_WELL in LIST ){ 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{ 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%2CwcIdentifier%2CdataMeasDate%2CdataValue&measUnits=ACFT&wcIdentifier=*U:1*&wdid=',C_WELL,'&apiKey=PwevUQJCStcYZfqrOYbuyztmNPlUJWby'),skip=2) %>% select(wdid,year=datameasuate,pumping=dataValue,wcIdentifier))} 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)) try(C_DAT <- read_csv(paste0('https://dwr.state.co.us/Rest/GET/api/v2/structures/divrec/divrecyear/?format=csv&dateFormat=dateOnly&fields=wdid%2CwcIdentifier%2CdataMeasDate%2CdataValue&measUnits=ACFT&wcIdentifier=*U:1*&wdid=',C_WELL),skip=2) %>% select(wdid,year=dataMeasDate,pumping=dataValue,wcIdentifier))
} }
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")}} 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() #closeAllConnections()

41
temp.r
View File

@ -3,27 +3,54 @@ library(fixest)
library("geosphere") library("geosphere")
DF <- readRDS("Data/Output_Data/Full_Data_Set.rds") %>% mutate(year=as.numeric(year)) %>% group_by(wdid) %>% mutate(SBD1=max(SBD1)) %>% ungroup DF <- readRDS("Data/Output_Data/Full_Data_Set.rds") %>% mutate(year=as.numeric(year)) %>% group_by(wdid) %>% mutate(SBD1=max(SBD1)) %>% ungroup
EMPTY_START <- DF %>% filter(year<=2010) %>% group_by(wdid) %>% filter(sum(AF)==0) %>% select(wdid) %>% unique TEMP <- read_csv("Data/Output_Data/Div3_Pumping_Data.csv")
TEMP %>% group_by(wdid,year) %>% filter(n()>1) %>% print(n=100)
#EMPTY_START <- DF %>% filter(year<=2010) %>% group_by(wdid) %>% filter(sum(AF)==0) %>% select(wdid) %>% unique
#DF <- DF %>% anti_join(EMPTY_START) #DF <- DF %>% anti_join(EMPTY_START)
DF[DF$CREP_any==1,'SBD1'] <- 1 OLD <- readRDS("Test/df_half.Rds") %>% mutate(wdid=GW_wdid,SBD1_OLD=SBD1,year=Year,AF_OLD=AF) %>% select(wdid,year,SBD1_OLD,AF_OLD)
DF[DF$CREP_any==1 &DF$SBD1_year==Inf,'SBD1_year'] <- 2011 BOTH <- DF %>% select(wdid,year,SBD1,AF) %>% full_join(OLD) %>% filter(year<2020)
BOTH %>% filter(is.na(SBD1))
BOTH %>% filter(AF!=AF_OLD) %>% print(n=800)
BOTH %>% filter(AF!=AF_OLD)
OLD %>% filter(wdid=='2405396')
BOTH %>% filter(wdid=='2405396')
BOTH %>% filter(is.na(SBD1_OLD))
DF$SBD1 <- ifelse(DF$SBD1_year!=Inf,1,0) DF$SBD1 <- ifelse(DF$SBD1_year!=Inf,1,0)
#DF <- DF %>% select(-ditch_2200627,-ditch_2200541,-ditch_3500570) #DF <- DF %>% select(-ditch_2200627,-ditch_2200541,-ditch_3500570)
DF$SBD1_year <- ifelse(DF$SBD1_year==2009,2011,DF$SBD1_year) DF$SBD1_year <- ifelse(DF$SBD1_year==2009,2011,DF$SBD1_year)
DF$POST <- ifelse(DF$year>=DF$SBD1_year,1,0) DF$POST <- ifelse(DF$year>=DF$SBD1_year,1,0)
DF <- DF %>% left_join(DF %>% group_by(wdid,CREP_any) %>% summarize(CREP_year=min(year)) %>% filter(CREP_any==1) %>% select(-CREP_any)) %>% mutate(CREP_year=ifelse(is.na(CREP_year),Inf,CREP_year)) DF <- DF %>% left_join(DF %>% group_by(wdid,CREP_any) %>% summarize(CREP_year=min(year)) %>% filter(CREP_any==1) %>% select(-CREP_any)) %>% mutate(CREP_year=ifelse(is.na(CREP_year),Inf,CREP_year))
AUG <- read_csv("Data/Augmentation/Augmentation.csv") %>% mutate(wdid=as.character(wdid)) %>% select(wdid)
DF <- DF %>% anti_join(AUG)
#feols(AF~SBD1*POST+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD2+SBD3+SBD4+SBD5+SBD6+wdid+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF) #feols(AF~SBD1*POST+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD2+SBD3+SBD4+SBD5+SBD6+wdid+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF)
########### ###########
#DF$SBD1_year <- ifelse(DF$SBD1_year==Inf,10000,DF$SBD1_year) #DF$SBD1_year <- ifelse(DF$SBD1_year==Inf,10000,DF$SBD1_year)
DF$POST_TEST <- ifelse(DF$year>=2011,1,0)
DF %>% filter(year==2009) %>% pull(wdid) %>% unique %>% length
TEST_DATA <- DF %>% mutate(SBD1=ifelse(SBD1==1 & SBD1_year<2020,1,0))
SUN_MOD <- feols(AF~SBD1+sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_alfalfa+per_potatoes|wdid+year+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,data=DF %>% filter(SBD1_year==Inf|SBD1_year==2011) ) SUN_MOD <- feols(AF~SBD1+sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_alfalfa+per_potatoes|wdid+year+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,data=DF %>% filter(SBD1_year==Inf|SBD1_year==2011) )
#coefplot(SUN_MOD) etable(feols(AF~SBD1*POST_TEST|wdid+year,data=TEST_DATA %>% filter(year<2019)))
coefplot(SUN_MOD)
SUN_MOD <- feols(AF~SBD1+sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028+per_alfalfa+per_potatoes|wdid+year+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,data=DF )
TEST <- DF %>% group_by(year,contacts,SBD1) %>% summarize(AF=sum(AF),SBD1_year=min(SBD1_year),POST=max(POST),CREP_temp=max(CREP_temp),CREP_perm=max(CREP_perm),SBD1_temp_fallow=max(SBD1_temp_fallow),SB22_028=max(SB22_028),SBD_half_fallow=max(SBD_half_fallow),SBD1_purchase=max(SBD1_purchase),SBD2=max(SBD2),SBD3=max(SBD3),SBD4=max(SBD4),SBD5=max(SBD5),SBD6=max(SBD6),ditch_2000631=max(ditch_2000631),ditch_Other=max(ditch_Other),ditch_2000829=max(ditch_2000829),ditch_2000798=max(ditch_2000798),ditch_2000816=max(ditch_2000816),ditch_2000753=max(ditch_2000753),ditch_2000753=max(ditch_2000753),ditch_2000623=max(ditch_2000623),CREP_any=max(CREP_any)) %>% ungroup
TEST %>% arrange(contacts,year)
TEST
TEST %>% select(contacts,year,SBD1,POST) %>% arrange(contacts,year) %>% filter(SBD1==1)
TEST_REG <- feols(log(AF+0.001)~sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|contacts+SBD2^year+SBD3^year+SBD4^year+SBD5^year+SBD6^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year,TEST )
TEST_REG <- feols(AF~sunab(SBD1_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|contacts,TEST )
coefplot(TEST_REG)
TEST
etable(TEST_REG ,agg="cohort")
coefplot(SUN_MOD)
#coefplot(feols(AF~sunab(CREP_year,year)|wdid+year,data=DF)) #coefplot(feols(AF~sunab(CREP_year,year)|wdid+year,data=DF))
coefplot(feols(AF~sunab(CREP_year,year)|wdid+year+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF)) coefplot(feols(AF~sunab(CREP_year,year)|wdid+year+ditch_2000812^year+ditch_2000631^year+ditch_Other^year+ditch_2000829^year+ditch_2000798^year+ditch_2000816^year+ditch_2000753^year+ditch_2000753^year+ditch_2000623^year+ditch_2200627^year+ditch_2200541^year+ditch_3500570^year,data=DF))