notebook.community

Edit and run

The optimize-waterallocation.jl script optimizes the joint extraction of surface and groundwater withdrawals and reservoir storage to satisfy given water demands. It produces four files:

withdrawals.jld: The amount withdrawn from each gauge along each canal.
returns.jld: The amount returned to each gauge along each canal.
captures.jld: The amount stored in each reservoir from one period to the next.
waterfromgw.jld: The amount of pumped water from the groundwater for each county.

How does the optimal division between ground and surface water differ from the recorded values?

First, let's take the recorded values.



In [2]:

    
using DataFrames
df = readtable("../data/extraction/USGS-2010.csv")[:, [:FIPS, :TO_SW, :TO_GW, :TO_To]]









    Out[2]:




FIPS TO_SW TO_GW TO_To
1 1001 37.17 18.38 55.55
2 1003 7.01 57.35 64.36
3 1005 4.34 6.71 11.05
4 1007 1.55 6.06 7.61
5 1009 53.04 4.04 57.08
6 1011 0.43 2.76 3.19
7 1013 0.79 3.73 4.52
8 1015 6.8 22.51 29.31
9 1017 4.45 0.95 5.4
10 1019 2.37 3.17 5.54
11 1021 1.97 5.35 7.32
12 1023 40.87 2.02 42.89
13 1025 21.25 3.06 24.31
14 1027 1.9 0.61 2.51
15 1029 0.88 0.91 1.79
16 1031 3.69 10.97 14.66
17 1033 1342.25 2.07 1344.32
18 1035 0.2 2.2 2.4
19 1037 0.09 0.84 0.93
20 1039 3.3 7.02 10.32
21 1041 0.64 2.52 3.16
22 1043 34.6 3.28 37.88
23 1045 1.61 8.64 10.25
24 1047 36.68 12.44 49.12
25 1049 7.9 4.04 11.94
26 1051 11.08 4.94 16.02
27 1053 34.59 10.62 45.21
28 1055 140.87 5.43 146.3
29 1057 3.99 0.66 4.65
30 1059 4.21 2.22 6.43
&vellip &vellip &vellip &vellip &vellip



In [3]:

    
surfaces = deserialize(open("../data/cache/counties/extraction/withdrawals.jld"));
grounds = deserialize(open("../data/cache/counties/extraction/waterfromgw.jld"))









    Out[3]:





3109x2 Array{Float64,2}:
 38401.5         38401.5      
 44491.8          3113.7      
  7638.82         7638.82     
  1230.13         5260.76     
     0.0         39459.2      
  2205.23         2205.23     
  3124.66            0.0      
 20202.1         20261.9      
  3733.0          3733.0      
  3829.78         3175.47     
  5060.29          442.333    
     0.0         29649.7      
 16805.4         16805.4      
     ⋮                        
     0.0             0.0      
     0.0         60323.1      
     0.0             0.0      
     0.0             0.0      
     0.0             0.0      
     0.0             1.16971e5
 72937.4        114681.0      
     1.12241e5   51872.9      
     0.0             0.0      
     0.0             0.0      
     0.0             0.0      
 19272.3         19280.2

Sum these to the county-year level.



In [12]:

    
drawsdata = read_rda("../data/countydraws.RData", convertdataframes=true);
draws = drawsdata["draws"]

optsw = []
optgw = []
for ii in 1:size(grounds, 1)
    push!(optsw, sum(surfaces[draws[:fips] .== df[ii, :FIPS], :]))
    push!(optgw, sum(grounds[ii, :]))
end

df[:optsw] = optsw / (1382592. / 1000.)
df[:optgw] = optgw / (1382592. / 1000.)
df









    Out[12]:




FIPS TO_SW TO_GW TO_To optsw optgw
1 1001 37.17 18.38 55.55 0.0 55.55
2 1003 7.01 57.35 64.36 29.927926355402754 34.43207364459724
3 1005 4.34 6.71 11.05 0.0 11.049999999999999
4 1007 1.55 6.06 7.61 2.9152733553631127 4.694726644636887
5 1009 53.04 4.04 57.08 28.54 28.54
6 1011 0.43 2.76 3.19 0.0 3.1899999999999995
7 1013 0.79 3.73 4.52 2.26 2.26
8 1015 6.8 22.51 29.31 0.0432667518520312 29.266733248147965
9 1017 4.45 0.95 5.4 0.0 5.4
10 1019 2.37 3.17 5.54 0.47324811286587337 5.066751887134127
11 1021 1.97 5.35 7.32 3.340069881408389 3.9799301185916103
12 1023 40.87 2.02 42.89 21.445 21.445
13 1025 21.25 3.06 24.31 0.0 24.309999999999995
14 1027 1.9 0.61 2.51 0.0 2.51
15 1029 0.88 0.91 1.79 0.0 1.7899999999999998
16 1031 3.69 10.97 14.66 10.407596469214596 4.252403530785405
17 1033 1342.25 2.07 1344.32 1293.2154137198188 51.1045862801811
18 1035 0.2 2.2 2.4 0.0 2.3999999999999995
19 1037 0.09 0.84 0.93 0.0 0.93
20 1039 3.3 7.02 10.32 5.159999999999999 5.159999999999999
21 1041 0.64 2.52 3.16 2.741862476481112 0.4181375235188876
22 1043 34.6 3.28 37.88 37.88 0.0
23 1045 1.61 8.64 10.25 0.8014402842357939 9.448559715764205
24 1047 36.68 12.44 49.12 24.818572263749846 24.301427736250144
25 1049 7.9 4.04 11.94 5.033198717937003 6.906801282062998
26 1051 11.08 4.94 16.02 0.17125394552160814 15.84874605447839
27 1053 34.59 10.62 45.21 4.577774980748596 40.6322250192514
28 1055 140.87 5.43 146.3 141.99449139367218 4.305508606327837
29 1057 3.99 0.66 4.65 0.0 4.65
30 1059 4.21 2.22 6.43 0.33271191415517626 6.0972880858448235
&vellip &vellip &vellip &vellip &vellip &vellip &vellip



In [13]:

    
include("../src/lib/graphing.jl")









    Out[13]:





usmap (generic function with 2 methods)

Where do we have under-allocations-- that is, where the optimized values do not reach the recorded levels?



In [15]:

    
df[:fips] = df[:FIPS]
df[:value] = convert(DataVector{Float64}, df[:TO_To] - df[:optsw] - df[:optgw])
usmap(df, true)









    












    Out[15]:





RCall.RObject{RCall.VecSxp}

Where do we use more or less GW than recorded (since this is the current source of prices)?



In [18]:

    
df[:value] = convert(DataVector{Float64}, (df[:TO_GW] - df[:optgw]))
usmap(df, true)









    












    Out[18]:





RCall.RObject{RCall.VecSxp}



In [ ]:

	FIPS	TO_SW	TO_GW	TO_To
1	1001	37.17	18.38	55.55
2	1003	7.01	57.35	64.36
3	1005	4.34	6.71	11.05
4	1007	1.55	6.06	7.61
5	1009	53.04	4.04	57.08
6	1011	0.43	2.76	3.19
7	1013	0.79	3.73	4.52
8	1015	6.8	22.51	29.31
9	1017	4.45	0.95	5.4
10	1019	2.37	3.17	5.54
11	1021	1.97	5.35	7.32
12	1023	40.87	2.02	42.89
13	1025	21.25	3.06	24.31
14	1027	1.9	0.61	2.51
15	1029	0.88	0.91	1.79
16	1031	3.69	10.97	14.66
17	1033	1342.25	2.07	1344.32
18	1035	0.2	2.2	2.4
19	1037	0.09	0.84	0.93
20	1039	3.3	7.02	10.32
21	1041	0.64	2.52	3.16
22	1043	34.6	3.28	37.88
23	1045	1.61	8.64	10.25
24	1047	36.68	12.44	49.12
25	1049	7.9	4.04	11.94
26	1051	11.08	4.94	16.02
27	1053	34.59	10.62	45.21
28	1055	140.87	5.43	146.3
29	1057	3.99	0.66	4.65
30	1059	4.21	2.22	6.43
&vellip	&vellip	&vellip	&vellip	&vellip