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

Updating file locations

Browse Source
This commit is contained in:
Alex Gebben Work 2026-07-22 16:25:13 -06:00
parent 8f386f0057
commit e2e4a32e0c
16 changed files with 167000 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.gitignore vendored
View File

@ -1,4 +1,13 @@
# ---> R
#Crop choice has very large raw files. These will be downloaded with a script so ignore them in git
Data/Crop_Choice/
#
Data/CREP_and_Fallow/Bill_SB22_Fallow_Payments.csv
Data/CREP_and_Fallow/CREP.csv
Data/CREP_and_Fallow/SBD1_Half_Fallow_Program.csv
Data/CREP_and_Fallow/SBD1_Temporary_Fallow_Payments.csv
Data/CREP_and_Fallow/SBD1_Well_Purchase_Program.csv
Data/Crop_Parcel_Data/Well_Ditch_Link.rds
*.swp
Data/ARP/*
Data/Input_Data/*

4
1_Create_CREP_Data.r
View File

@ -122,6 +122,7 @@ CURRENT <- IN_PROGRAM(i,SBD1_TEMP_FALLOW)
if(exists("ALL_PROGRAM_WELL_DATA")){ALL_PROGRAM_WELL_DATA <- rbind(ALL_PROGRAM_WELL_DATA,CURRENT)} else{ALL_PROGRAM_WELL_DATA <- CURRENT}
}
PROGRAM_ALL <- function(DATA_SET){do.call(rbind,lapply(DATA_SET$wdid,function(x){IN_PROGRAM(x,DATA_SET)})) }
ALL_PROGRAMS <- rbind(PROGRAM_ALL(CREP_PERM),
PROGRAM_ALL(CREP_TEMP) %>% mutate(program='CREP_temp'),
@ -134,5 +135,4 @@ ALL_PROGRAMS <- ALL_PROGRAMS %>% mutate(CREP_any=ifelse(CREP_perm+CREP_temp>0,1,
dir.create("Data/Output_Data",showWarnings=FALSE,recursive=TRUE)
write_csv(ALL_PROGRAMS,"Data/Output_Data/Fallow_Program_Data.csv")
saveRDS(ALL_PROGRAMS,"Data/Output_Data/Fallow_Program_Data.rds")
print("Script 1: Create CREP data completed")

41
Data/CREP_and_Fallow/Create_CREP_Data.r → 2_Create_Fallow_Program_Data.r
View File

@ -1,17 +1,26 @@
library(tidyverse)
library(janitor)
#Data manually copied from tables in the yearly reports, and combined here. The original PDF files are very large but can be reviewed for accuracy by going to the public P-drive link to the data.
#The link to this raw data is https://u.pcloud.link/publink/show?code=kZ0vRc5Zv5LhTm5PoBbj5lyoSDgwwfX0DUCy or it can be found in the Annual Reports on the Subdistrict website.
#As an example the file "ARP 2023- ALL for Website.pdf" in the drive has a table in Appendix L- Crep and Fallow Programs starting at page 314. Each column was manually copied into the vectors in this file before being converted into a data frame and saved for later use in the data collation, and final analysis.
###########Some data was copied into rough csv files from the pdf's to make data maniputlation easier. This is the location of those files
ROOT_DIR <- "./Data/CREP_and_Fallow/"
dir.create(ROOT_DIR,showWarnings=FALSE)
RAW_DATA_DIR <- paste0(ROOT_DIR,"Raw_CSV_Data_Made_from_Files/")
################Perm CREP contract
CONTRACT_NUMBER <- c('ALA#3','ALA#6','ALA#7','ALA#8','ALA#9','ALA#10','ALA#12','ALA#15','SAG#6','ALA#17','ALA#18','ALA#22','ALA#23','ALA#25','RG#4','ALA#26','ALA#27','ALA#28','ALA#29','ALA#30','ALA#31','ALA#32','ALA#33','ALA#34','ALA#38','ALA#39','SAG#33','SAG#34','ALA#40','ALA#41','ALA#42','ALA#43','ALA#44','ALA#45','ALA#47','SAG#39')
FIRST_FALLOW_YEAR <- c(2014,2014,2014,2014,2014,2014,2014,2014,2015,2015,2015,2015,2015,2015,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2018,2018,2020,2020,2020,2020,2020,2021,2021,2021,2021,2024)
ACRES <- c(124.9,126,119.5,119.2,121.1,118.1,122.8,67,114.1,118.6,122,121,124.66,80,149.8,110,110,110,92.9,122.3,94,123,126,126,121.28,120.5,122.8,122,118,121.84,120,120,120.01,120.1,120.11,122.94)
CREP_PERM <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,ACRES,read_csv("CREP_Perm.csv")) %>% as_tibble %>% mutate(program='CREP',RETURN_YEAR=Inf,contract_type='Perm')
CREP_PERM <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,ACRES,read_csv(paste0(RAW_DATA_DIR,"CREP_Perm.csv"))) %>% as_tibble %>% mutate(program='CREP',RETURN_YEAR=Inf,contract_type='Perm')
###############################Temp CREP contract
CONTRACT_NUMBER <- c('SAG#1','SAG#2','SAG#3','SAG#4','RG#1','RG#2','ALA#2','ALA#11','SAG#7','SAG#8','SAG#9','SAG#10','SAG#11','ALA#16','ALA#19','ALA#21','ALA#24','SAG#12','SAG#13','RG#3','RG#7','RG#8','ALA#35','SAG#14','SAG#15','SAG#16','ALA#36','SAG#17','SAG#18','SAG#19','SAG#20','SAG#21','SAG#22','SAG#23','SAG#24','SAG#25','SAG#26','SAG#27','SAG#28','ALA#37','SAG#29','SAG#30','SAG#31','RG#9','RG#10','SAG#32','RG#11','SAG#35','SAG#36','SAG#37','SAG#38','RG#12','RG#13')
FIRST_FALLOW_YEAR <- c(2014,2014,2014,2014,2014,2014,2014,2014,2015,2015,2015,2015,2015,2015,2015,2015,2015,2016,2016,2016,2016,2016,2016,2017,2017,2017,2017,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2018,2019,2019,2019,2019,2019,2020,2022,2023,2023,2023,2023,2024,2024)
RETURN_YEAR <- c(2029,2029,2029,2029,2029,2029,2029,2029,2030,2030,2030,2030,2030,2030,2030,2030,2030,2031,2031,2031,2031,2031,2031,2032,2032,2032,2032,2033,2033,2033,2033,2033,2033,2033,2033,2033,2033,2033,2033,2033,2034,2034,2034,2034,2034,2035,2037,2038,2038,2038,2038,2038,2038)
ACRES <- c(144,144,210,60,130,120.4,120,121.5,172.09,113,191,116.5,120,124,120,129,120.97,120,124,139.9,122,123.32,122,120,122.4,123.4,113.92,120,120.35,114.32,124.78,125.58,119.3,123,125.15,126.1,126.3,125.5,53.6,106,112.81,126.95,118.9,118.36,120,120,100,120,130.02,114.54,120,121.92,126.15)
CREP_TEMP <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,read_csv("CREP_Temp.csv")) %>% as_tibble %>% mutate(program='CREP',contract_type='Temp')
CREP_TEMP <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,read_csv(paste0(RAW_DATA_DIR,"CREP_Temp.csv"))) %>% as_tibble %>% mutate(program='CREP',contract_type='Temp')
######SBD1 Fallow Program
CONTRACT_NUMBER <- c(paste("Fallow Parcel ",1:23),paste("Fallow Parcel ",1:27),c(paste("Fallow Parcel ",1:9),paste("Fallow Parcel ",11:33)))
@ -40,7 +49,7 @@ ACRES <- c(
)
SBD1_TEMP_FALLOW <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,read_csv("SBD1_Temp_Fallow.csv")) %>% as_tibble %>% mutate(program='SBD Temp Fallow',contract_type='Temp',wdid1=as.character(wdid1),wdid2=as.character(wdid2),wdid3=as.character(wdid3),wdid4=as.character(wdid4),wdid5=as.character(wdid5),wdid6=as.character(wdid6),wdid7=as.character(wdid7),wdid8=as.character(wdid8))
SBD1_TEMP_FALLOW <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,read_csv(paste0(RAW_DATA_DIR,"SBD1_Temp_Fallow.csv"))) %>% as_tibble %>% mutate(program='SBD Temp Fallow',contract_type='Temp',wdid1=as.character(wdid1),wdid2=as.character(wdid2),wdid3=as.character(wdid3),wdid4=as.character(wdid4),wdid5=as.character(wdid5),wdid6=as.character(wdid6),wdid7=as.character(wdid7),wdid8=as.character(wdid8))
SBD1_TEMP_FALLOW
nrow(SBD1_TEMP_FALLOW )
(SBD1_TEMP_FALLOW %>% filter(ACRES==125))[,-1:-1]
@ -65,19 +74,18 @@ FIRST_FALLOW_YEAR <- c(rep(2021,4),rep(2024,42))
RETURN_YEAR <- rep(2025,46)
ACRES <- c(73.46,50,119,118,122.01,126,126,120,198,122,123,116,126,126,120,64,124,116,126,126,126,126,120,120,126,126,100,126,126,118.62,139.62,101.4,120,118.8,118.32,118,117,136,120,114.68,114.67,49,126.74,120,126,120)
DATA_2024 <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,read_csv("SBD1_Temp_Fallow_2024.csv")) %>% as_tibble %>% mutate(program='SBD Temp Fallow',contract_type='Temp',wdid1=as.character(wdid1),wdid2=as.character(wdid2),wdid3=as.character(wdid3),wdid4=as.character(wdid4),wdid5=as.character(wdid5),wdid6=as.character(wdid6))
DATA_2024 <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,read_csv(paste0(RAW_DATA_DIR,"SBD1_Temp_Fallow_2024.csv"))) %>% as_tibble %>% mutate(program='SBD Temp Fallow',contract_type='Temp',wdid1=as.character(wdid1),wdid2=as.character(wdid2),wdid3=as.character(wdid3),wdid4=as.character(wdid4),wdid5=as.character(wdid5),wdid6=as.character(wdid6))
SBD1_TEMP_FALLOW <- SBD1_TEMP_FALLOW %>% full_join(DATA_2024)
####SBD1 Permanent well Retirement
CONTRACT_NUMBER <- c('2021-01','2021-02','2021-03','2021-04','2021-05','2021-06','2021-07','2021-08','2021-09','2021-10','2021-11','2022-1-01','2022-1-02','2022-1-03','2022-1-04','2022-1-05','2022-1-06','2022-1-07','2022-1-08','2022-2-01','2023-01','2023-02','2023-03','2023-05','2023-06','2023-07','2023-08','2023-09')
FIRST_FALLOW_YEAR <- c(2022,2022,2022,2022,2022,2022,2022,2022,2022,2022,2022,2023,2023,2023,2023,2023,2023,2023,2023,2023,2024,2024,2024,2024,2024,2024,2024,2024)
ACRES <- c(125,135,123,130,117,122,97,121,125,125,128,130,128,119,124,121,124,121,124,123,124,130,129,123,120,240,123,118)
SBD1_WELL_PURCHASE <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,ACRES,read_csv("SBD1_Well_Purchase.csv")) %>% as_tibble %>% mutate(program='SBD1 Purchase',RETURN_YEAR=Inf,contract_type='Perm',Data_Year=2025)
SBD1_WELL_PURCHASE
SBD1_WELL_PURCHASE <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,ACRES,read_csv(paste0(RAW_DATA_DIR,"SBD1_Well_Purchase.csv"))) %>% as_tibble %>% mutate(program='SBD1 Purchase',RETURN_YEAR=Inf,contract_type='Perm',Data_Year=2025)
######Groundwater Compact Compliance and Sustainability Fund (SB22-028)
CONTRACT_NUMBER <- c('004','008','009','009','011','011','117','118','119','120','121','121','121','121','122','123','125','126','231','231','231','233','235','239','340','340','342')
FIRST_FALLOW_YEAR <- c(2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2024,2025,2025,2025,2025)
ACRES <- c(125,118,121,131,119,119,125,125,123,122,120,130,123,126,125,120,124,117,126,120,125,122,124,120,119,125,122)
SB22_Data <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,ACRES,read_csv("SB22_Data.csv")) %>% as_tibble %>% mutate(program='SB22-028',RETURN_YEAR=Inf,contract_type='Perm',Data_Year=2025)
SB22_Data <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,ACRES,read_csv(paste0(RAW_DATA_DIR,"SB22_Data.csv"))) %>% as_tibble %>% mutate(program='SB22-028',RETURN_YEAR=Inf,contract_type='Perm',Data_Year=2025)
###########Half Fallow Program
CONTRACT_NUMBER <- paste("Half Usage",1:38)
@ -85,7 +93,7 @@ FIRST_FALLOW_YEAR <- rep(2020,38)
RETURN_YEAR <- rep(2021,38)
LIMIT <- c(90.6,91.5,83.4,85.2,89.8,82.4,81.0,72.0,98.8,97.8,92.2,94.0,111.7,118.0,102.0,92.1,117.0,99.4,95.7,149.0,54.6,22.6,102.4,97.5,211.6,104.5,47.5,76.2,34.5,74.1,53.0,85.7,130.7,97.4,115.0,56.0,33.4,45.9)
ACRES <- c(111.92,119.22,116.82,117.26,119.1,123,118.03,106.98,109.15,120,117,128,126,126,125,121,121,126,121,120,43,24,138,137,117,121.5,142,126,174.22,118,118,119,126,126,120,114,126,63)
HALF_FALLOW <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,LIMIT,read_csv("SBD1_Half_Fallow.csv")) %>% as_tibble %>% mutate(program='SBD Half Fallow',contract_type='Temp')
HALF_FALLOW <- cbind(CONTRACT_NUMBER,FIRST_FALLOW_YEAR,RETURN_YEAR,ACRES,LIMIT,read_csv(paste0(RAW_DATA_DIR,"SBD1_Half_Fallow.csv"))) %>% as_tibble %>% mutate(program='SBD Half Fallow',contract_type='Temp')
#######################Combine
SBD1_TEMP_FALLOW$CONTRACT_NUMBER <- paste0(SBD1_TEMP_FALLOW$Data_Year,"_",SBD1_TEMP_FALLOW$CONTRACT_NUMBER)
SBD1_TEMP_FALLOW <- SBD1_TEMP_FALLOW %>% pivot_longer(c(wdid1,wdid2,wdid3,wdid4,wdid5,wdid6,wdid7,wdid8),values_to='wdid') %>% select(-name) %>% filter(!is.na(wdid))%>% clean_names
@ -124,18 +132,15 @@ if(exists("ALL_PROGRAM_WELL_DATA")){ALL_PROGRAM_WELL_DATA <- rbind(ALL_PROGRAM
}
ALL_PROGRAM_WELL_DATA %>% print(n=100)
IN_PROGRAM(2013956,CREP_PERM)
IN_PROGRAM(2705126,CREP_TEMP)
#IN_PROGRAM(2013956,CREP_PERM)
#IN_PROGRAM(2705126,CREP_TEMP)
dir.create("./Cleaned",showWarnings=FALSE)
write_csv(SBD1_TEMP_FALLOW,"./Cleaned/SBD1_Temporary_Fallow_Payments.csv")
write_csv(SBD1_WELL_PURCHASE,"./Cleaned/SBD1_Well_Purchase_Program.csv")
write_csv(HALF_FALLOW,"./Cleaned/SBD1_Half_Fallow_Program.csv")
write_csv(SB22_Data,"./Cleaned/Bill_SB22_Fallow_Payments.csv")
write_csv(CREP,"./Cleaned/CREP.csv")
write_csv(SBD1_TEMP_FALLOW,paste0(ROOT_DIR,"SBD1_Temporary_Fallow_Payments.csv"))
write_csv(SBD1_WELL_PURCHASE,paste0(ROOT_DIR,"SBD1_Well_Purchase_Program.csv"))
write_csv(HALF_FALLOW,paste0(ROOT_DIR,"SBD1_Half_Fallow_Program.csv"))
write_csv(SB22_Data,paste0(ROOT_DIR,"Bill_SB22_Fallow_Payments.csv"))
write_csv(CREP,paste0(ROOT_DIR,"CREP.csv"))

0
3_Proc_Dynamic_Data.r → 3_Download_and_Process_Pumping_Data.r
View File

6
Data/Crop_Parcel_Data/Proc_Ditches.r → 3_Process_Ditches.r
View File

@ -1,5 +1,5 @@
library(tidyverse)
ALL_PARCEL_DATA <- read_csv("WELL_DITCH_PARCEL.csv")
ALL_PARCEL_DATA <- read_csv("Data/Crop_Parcel_Data/WELL_DITCH_PARCEL.csv")
NAMES <- colnames(ALL_PARCEL_DATA )
PARCEL_LINK <- ALL_PARCEL_DATA %>% select(NAMES[grep("PARCEL_ID|GW_ID|SW_WDID",NAMES )])
PARCEL_LINK <- PARCEL_LINK %>% pivot_longer(-PARCEL_ID,values_to='wdid',names_to="type") %>% filter(!is.na(wdid))
@ -17,4 +17,6 @@ DITCH_WELL <- DITCH_WELL %>% pivot_wider(values_from=ditch_id,names_from=ditch_i
DITCH_WELL[,-1] <- ifelse(is.na(DITCH_WELL[,-1]),0,1)
DITCH_WELL <- DITCH_WELL %>% mutate(wdid=as.character(wdid))
saveRDS(DITCH_WELL,"Well_Ditch_Link.rds")
saveRDS(DITCH_WELL,"Data/Crop_Parcel_Data/Well_Ditch_Link.rds")
print("Script 2: Process ditch link data completed")

12
4_Download_Crop_Choice_Data_frmo_Cloud.r Normal file
View File

@ -0,0 +1,12 @@
library(RCurl)
#Download very large crop choice files created using the Hyrdobase map file in QGIS. This has every crop and technology combination for all parcels in 2002, or 2005. This is used as the pre-treatment control for crop choice correlated with factors such as soil quality.
#These files are excluded from git due to the very large space requirment. Instead Alex Gebben has hosted them on a Pcloud drive, and provided public access, allowing git to ignore them but have R download at project start.
#Location of the files
CROP_2002_URL <- 'https://def3.pcloud.com/DLZHyJ74J7ZfHere67ZCPjOZXZqxJG5kZ2ZZM3FZZTvmJZzYZsLZ3TZH807r8lSpP0MjbM4PoY9z7sxqCDy/IRRIG_2002.csv'
CROP_2005_URL <- 'https://def3.pcloud.com/DLZwyJ74J7Z7zere67ZCPjOZXZExJG5kZ2ZZM3FZZuVFHZ3YZnYZmgZ9u59W63pgNRX6L94jxED10e4TrYk/IRRIG_2005.csv'
DEST_DIR <- "./Data/Crop_Choice/"
dir.create(DEST_DIR,showWarnings=FALSE)
download.file(CROP_2002_URL,destfile=paste0(DEST_DIR,"IRRIG_2002.csv"))
download.file(CROP_2005_URL,destfile=paste0(DEST_DIR,"IRRIG_2005.csv"))

3
2_Proc.r → 5_Process_All_Cleaned_Data.r
View File

@ -38,7 +38,7 @@ STATIC_DATA[,DITCH_COL_NAMES] <- STATIC_DATA[,DITCH_COL_NAMES]%>% replace(is.na(
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 no crop data was available.
STATIC_DATA <- STATIC_DATA %>% replace(is.na(.), 0)
STATIC_DATA[-1:-5] <- STATIC_DATA[-1:-5] %>% replace(is.na(.), 0)
###Make sure all CREP wells are included in SBD1
PROGRAM_WELLS <- c(read_csv("Data/CREP_and_Fallow/Cleaned/Fallow_Program_Data/Bill_SB22_Fallow_Payments.csv")$wdid,
read_csv("Data/CREP_and_Fallow/Cleaned/Fallow_Program_Data/CREP.csv")$wdid,
@ -80,4 +80,5 @@ PUMPING <- PUMPING %>% pivot_wider(values_from=AF,names_from=year)%>% group_by(w
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")
print("Script 3: Process all data into panel completed")

0
Data/CREP_and_Fallow/CREP_Perm.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/CREP_Perm.csv
View File

0
Data/CREP_and_Fallow/CREP_Temp.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/CREP_Temp.csv
View File

0
Data/CREP_and_Fallow/SB22_Data.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/SB22_Data.csv
View File

Can't render this file because it contains an unexpected character in line 3 and column 9.

0
Data/CREP_and_Fallow/SBD1_Half_Fallow.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/SBD1_Half_Fallow.csv
View File

Can't render this file because it has a wrong number of fields in line 2.

0
Data/CREP_and_Fallow/SBD1_Temp_Fallow.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/SBD1_Temp_Fallow.csv
View File

Can't render this file because it has a wrong number of fields in line 2.

0
Data/CREP_and_Fallow/SBD1_Temp_Fallow_2024.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/SBD1_Temp_Fallow_2024.csv
View File

Can't render this file because it has a wrong number of fields in line 2.

0
Data/CREP_and_Fallow/SBD1_Well_Purchase.csv → Data/CREP_and_Fallow/Raw_CSV_Data_Made_from_Files/SBD1_Well_Purchase.csv
View File

Can't render this file because it has a wrong number of fields in line 2.

166948
Data/Crop_Parcel_Data/WELL_DITCH_PARCEL.csv Normal file
View File

File diff suppressed because it is too large Load Diff

115
temp.r
View File

@ -1,115 +0,0 @@
library(tidyverse)
library(fixest)
library("geosphere")
DF <- readRDS("Data/Output_Data/Full_Data_Set.rds") %>% mutate(year=as.numeric(year)) %>% group_by(wdid) %>% mutate(SBD1=max(SBD1)) %>% ungroup
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)
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)
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 <- DF %>% select(-ditch_2200627,-ditch_2200541,-ditch_3500570)
DF$SBD1_year <- ifelse(DF$SBD1_year==2009,2011,DF$SBD1_year)
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))
#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$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) )
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+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))
GET_PROGRAM <-function(PROGRAM_NUMBER,DATA,CUTOFF=0.25,YEAR_RANGE=2009:2025){
#DATA <- DF
#PROGRAM_NUMBER <- 4
#YEAR_RANGE=2009:2025
#CUTOFF <- 0.25
DATA_ORIG <- DATA
MASTER_WELL_LIST <- DATA %>% select(wdid,longitude,latitude) %>% unique
WELLS_IN_PROGRAM <- DATA[DATA[,PROGRAM_NUMBER]==1,]
PROGRAM_START <- min(WELLS_IN_PROGRAM$year)
PROGRAM_YEAR_RANGE <- min(WELLS_IN_PROGRAM$year):max(YEAR_RANGE)
for(YEAR in PROGRAM_YEAR_RANGE){
# YEAR <-2014
C_PROGRAM_WELLS <- WELLS_IN_PROGRAM[WELLS_IN_PROGRAM[,"year"]==YEAR,c("longitude","latitude")] %>% unique %>% as.matrix
RES_DISTANCE <- c()
for(i in 1:nrow(MASTER_WELL_LIST)){
RES_DISTANCE[i] <-ifelse((2.38084*min(distHaversine(MASTER_WELL_LIST[i,c("longitude","latitude")],C_PROGRAM_WELLS))/5280)<=CUTOFF,1,0)
}
C_RES <- cbind(MASTER_WELL_LIST[,1],RES_DISTANCE) %>% as_tibble
colnames(C_RES) <-c("wdid",paste0("close_",colnames(DATA_ORIG)[PROGRAM_NUMBER]))
C_RES$year <- YEAR
if(exists("RES")){RES <- rbind(RES,C_RES)} else{RES <- C_RES}
}
if(PROGRAM_START>min(YEAR_RANGE)){
FILL_IN_YEARS <- min(YEAR_RANGE):(PROGRAM_START-1)
length(FILL_IN_YEARS )
5*(nrow(MASTER_WELL_LIST))
MASTER_WELL_LIST[,1]
FILL_IN_WELLS <- rep(t(MASTER_WELL_LIST[,1]),length(FILL_IN_YEARS))
ZEROS <- rep(0,length(FILL_IN_WELLS))
FILL_IN_YEARS <- rep( FILL_IN_YEARS,nrow(MASTER_WELL_LIST)) %>% sort
FILL_IN <- cbind(FILL_IN_WELLS,ZEROS,FILL_IN_YEARS) %>% as_tibble
FILL_IN <- FILL_IN %>% mutate(ZEROS=as.numeric(ZEROS),FILL_IN_YEARS=as.numeric(FILL_IN_YEARS))
colnames(FILL_IN) <- c("wdid",paste0("close_",colnames(DATA_ORIG)[PROGRAM_NUMBER]),"year")
RES <- rbind(FILL_IN,RES) %>% unique
}
RES <- RES[,c(1,3,2)]
return(RES)
}
DF <- DF %>% left_join(GET_PROGRAM(4,DF)) %>%left_join(GET_PROGRAM(5,DF)) %>% left_join(GET_PROGRAM(6,DF) )
CLOSE_YEAR <-
CLOSE_YEAR <-DF %>% filter(close_CREP_any==1) %>% group_by(wdid) %>% summarize(close_CREP_year=min(year)) %>% unique
DF <- DF %>% left_join(CLOSE_YEAR) %>% mutate(close_CREP_year=ifelse(is.na(close_CREP_year),Inf,close_CREP_year))
summary(DF$close_CREP_year)
coefplot(feols(AF~sunab(close_CREP_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD1^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 ))
etable(feols(AF~sunab(close_CREP_year,year)+CREP_temp+CREP_perm+SBD_half_fallow+SBD1_purchase+SBD1_temp_fallow+SB22_028|wdid+SBD1^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 ),agg="cohort")