optimize-waterallocation


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]:
FIPSTO_SWTO_GWTO_To
1100137.1718.3855.55
210037.0157.3564.36
310054.346.7111.05
410071.556.067.61
5100953.044.0457.08
610110.432.763.19
710130.793.734.52
810156.822.5129.31
910174.450.955.4
1010192.373.175.54
1110211.975.357.32
12102340.872.0242.89
13102521.253.0624.31
1410271.90.612.51
1510290.880.911.79
1610313.6910.9714.66
1710331342.252.071344.32
1810350.22.22.4
1910370.090.840.93
2010393.37.0210.32
2110410.642.523.16
22104334.63.2837.88
2310451.618.6410.25
24104736.6812.4449.12
2510497.94.0411.94
26105111.084.9416.02
27105334.5910.6245.21
281055140.875.43146.3
2910573.990.664.65
3010594.212.226.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]:
FIPSTO_SWTO_GWTO_Tooptswoptgw
1100137.1718.3855.550.055.55
210037.0157.3564.3629.92792635540275434.43207364459724
310054.346.7111.050.011.049999999999999
410071.556.067.612.91527335536311274.694726644636887
5100953.044.0457.0828.5428.54
610110.432.763.190.03.1899999999999995
710130.793.734.522.262.26
810156.822.5129.310.043266751852031229.266733248147965
910174.450.955.40.05.4
1010192.373.175.540.473248112865873375.066751887134127
1110211.975.357.323.3400698814083893.9799301185916103
12102340.872.0242.8921.44521.445
13102521.253.0624.310.024.309999999999995
1410271.90.612.510.02.51
1510290.880.911.790.01.7899999999999998
1610313.6910.9714.6610.4075964692145964.252403530785405
1710331342.252.071344.321293.215413719818851.1045862801811
1810350.22.22.40.02.3999999999999995
1910370.090.840.930.00.93
2010393.37.0210.325.1599999999999995.159999999999999
2110410.642.523.162.7418624764811120.4181375235188876
22104334.63.2837.8837.880.0
2310451.618.6410.250.80144028423579399.448559715764205
24104736.6812.4449.1224.81857226374984624.301427736250144
2510497.94.0411.945.0331987179370036.906801282062998
26105111.084.9416.020.1712539455216081415.84874605447839
27105334.5910.6245.214.57777498074859640.6322250192514
281055140.875.43146.3141.994491393672184.305508606327837
2910573.990.664.650.04.65
3010594.212.226.430.332711914155176266.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 [ ]: