compare-simulate


Here we compare the results of the simulation to recorded extractions from the USGS.

Total Irrigation Demand

Irrigation demand is estimated by the model as the withdrawals necessary to offset water deficits in irrigated land. That is, irrigated land is assumed to be completely and exactly irrigated.

First, we perform a complete simulation of the model.


In [1]:
include("../src/simulate.jl")


Loading from saved region network...
WARNING: imported binding for edges overwritten in module Main
Loading from saved water network...
Loading from saved region network...
Creating model...
  8.627338 seconds (2.32 M allocations: 134.071 MB, 1.11% gc time)

The results, stored in the model m, include irrigation demand, in units of $1000 m^3/year$.


In [2]:
simulated = getdataframe(m, :Agriculture, :totalirrigation)


Out[2]:
regionstimetotalirrigation
10100120150.0
20100320150.0
30100520150.0
40100720150.0
50100920150.0
60101120150.0
70101320150.0
80101520150.0
901017201513.337729972292596
100101920150.0
110102120150.0
120102320150.0
130102520150.0
140102720150.0
15010292015150.22831707787265
16010312015826.8036650819048
170103320150.0
180103520150.0
190103720150.0
200103920150.0
210104120151156.1243649795708
220104320150.0
230104520150.0
24010472015194.59543588021245
2501049201547.884290665963064
26010512015586.128567044208
2701053201515.413981383059623
280105520150.0
290105720150.0
300105920150.0
&vellip&vellip&vellip&vellip

Now we can compare these to recorded USGS irrigation. Start by loading the loading the USGS data, which is exported from an Excel file in prepare/extraction.


In [2]:
recorded = readtable("../data/extraction/USGS-2010.csv")


Out[2]:
STATESTATEFIPSCOUNTYCOUNTYFIPSFIPSYEARTP_TotPopPS_GWPS_SWPS_ToDO_GWDO_SWDO_ToIN_GWIN_SWIN_ToIR_GWIR_SWIR_ToLI_GWLI_SWLI_ToMI_GWMI_SWMI_ToPT_GWPT_SWPT_ToTO_GWTO_SWTO_To
1AL1Autauga County11001201054.5715.090.05.090.370.00.379.8931.1441.032.880.062.940.060.090.150.090.040.130.05.845.8418.3837.1755.55
2AL1Baldwin County310032010182.26522.970.022.971.710.01.710.00.00.032.36.6638.960.160.350.510.210.00.210.00.00.057.357.0164.36
3AL1Barbour County51005201027.4574.150.04.150.150.00.151.580.01.580.461.712.170.242.572.810.130.060.190.00.00.06.714.3411.05
4AL1Bibb County71007201022.9154.890.04.890.150.00.150.00.00.00.630.060.690.031.491.520.360.00.360.00.00.06.061.557.61
5AL1Blount County91009201057.3222.4452.1754.610.890.00.890.00.00.00.120.360.480.470.510.980.120.00.120.00.00.04.0453.0457.08
6AL1Bullock County111011201010.9142.30.02.30.080.00.080.00.00.00.290.320.610.050.090.140.040.020.060.00.00.02.760.433.19
7AL1Butler County131013201020.9472.70.02.70.280.00.280.30.00.30.020.320.340.430.470.90.00.00.00.00.00.03.730.794.52
8AL1Calhoun County1510152010118.57220.832.4723.30.510.00.510.970.00.970.04.094.090.130.210.340.070.030.10.00.00.022.516.829.31
9AL1Chambers County171017201034.2150.04.314.310.710.00.710.00.00.00.170.030.20.070.110.180.00.00.00.00.00.00.954.455.4
10AL1Cherokee County191019201025.9892.530.963.490.520.00.520.00.00.00.01.241.240.110.170.280.010.00.010.00.00.03.172.375.54
11AL1Chilton County211021201043.6433.061.834.890.740.00.740.350.00.351.120.031.150.070.110.180.010.00.010.00.00.05.351.977.32
12AL1Choctaw County231023201013.8591.360.01.360.620.00.620.040.7640.760.00.030.030.040.080.120.00.00.00.00.00.02.0240.8742.89
13AL1Clarke County251025201025.8332.250.93.150.380.00.380.020.2220.220.020.050.070.060.080.140.350.00.350.00.00.03.0621.2524.31
14AL1Clay County271027201013.9320.01.661.660.480.00.480.00.00.00.00.060.060.130.180.310.00.00.00.00.00.00.611.92.51
15AL1Cleburne County291029201014.9720.00.560.560.770.00.770.00.00.00.00.160.160.140.160.30.00.00.00.00.00.00.910.881.79
16AL1Coffee County311031201049.9487.60.07.60.720.00.721.130.01.130.72.973.670.820.721.540.00.00.00.00.00.010.973.6914.66
17AL1Colbert County331033201054.4280.578.228.790.270.00.270.2169.5469.750.891.32.190.130.160.290.00.730.730.01262.31262.32.071342.251344.32
18AL1Conecuh County351035201013.2281.690.01.690.320.00.320.00.00.00.060.030.090.130.170.30.00.00.00.00.00.02.20.22.4
19AL1Coosa County371037201011.5390.30.00.30.430.00.430.00.00.00.010.020.030.020.030.050.080.040.120.00.00.00.840.090.93
20AL1Covington County391039201037.7654.950.04.950.880.00.880.050.00.050.791.171.960.30.390.690.050.00.050.01.741.747.023.310.32
21AL1Crenshaw County411041201013.9062.050.02.050.190.00.190.00.00.00.00.250.250.280.390.670.00.00.00.00.00.02.520.643.16
22AL1Cullman County431043201080.4060.4830.5731.050.240.00.240.431.832.261.081.172.251.041.012.050.010.020.030.00.00.03.2834.637.88
23AL1Dale County451045201050.2516.850.06.850.690.00.690.00.00.00.871.292.160.180.290.470.050.030.080.00.00.08.641.6110.25
24AL1Dallas County471047201043.825.90.05.90.640.00.640.1432.1932.330.372.052.425.122.327.440.270.120.390.00.00.012.4436.6849.12
25AL1Dekalb County491049201071.1090.816.37.111.320.01.320.770.00.770.10.480.580.971.092.060.070.030.10.00.00.04.047.911.94
26AL1Elmore County511051201079.3033.79.6513.350.420.00.420.00.00.00.421.141.560.080.140.220.320.150.470.00.00.04.9411.0816.02
27AL1Escambia County531053201038.3195.680.05.680.630.00.631.4733.6635.132.180.843.020.080.090.170.580.00.580.00.00.010.6234.5945.21
28AL1Etowah County5510552010104.434.6915.6820.370.310.00.310.09.219.210.00.670.670.190.540.730.240.110.350.0114.66114.665.43140.87146.3
29AL1Fayette County571057201017.2410.051.982.030.520.00.520.01.11.10.00.170.170.090.110.20.00.630.630.00.00.00.663.994.65
30AL1Franklin County591059201031.7041.083.494.570.510.00.510.00.00.00.00.150.150.320.430.750.310.140.450.00.00.02.224.216.43
&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip&vellip

We will use R to display the result. First, let's map the irrigation as recorded by USGS. We need to prepare the shapefile first, by adding the FIPS code to each polygon point.


In [3]:
using RCall
R"library(ggplot2)"
R"library(PBSmapping)"
R"shapes <- importShapefile('../data/mapping/US_county_2000-simple')"
R"polydata <- attributes(shapes)$PolyData"
R"polydata$STATE <- as.numeric(levels(polydata$STATE))[polydata$STATE]"
R"polydata$COUNTY <- as.numeric(levels(polydata$COUNTY))[polydata$COUNTY]"
R"shapes$id <- polydata$STATE[shapes$PID] * 100 + polydata$COUNTY[shapes$PID] / 10";
R"names(shapes) <- tolower(names(shapes))";


WARNING: RCall.jl Warning: package ‘ggplot2’ was built under R version 3.2.3
WARNING: RCall.jl 
-----------------------------------------------------------
PBS Mapping 2.69.76 -- Copyright (C) 2003-2016 Fisheries and Oceans Canada

PBS Mapping comes with ABSOLUTELY NO WARRANTY;
for details see the file COPYING.
This is free software, and you are welcome to redistribute
it under certain conditions, as outlined in the above file.

A complete user guide 'PBSmapping-UG.pdf' is located at 
/Library/Frameworks/R.framework/Versions/3.2/Resources/library/PBSmapping/doc/PBSmapping-UG.pdf

Packaged on 2015-04-23
Pacific Biological Station, Nanaimo

All available PBS packages can be found at
http://code.google.com/p/pbs-software/

To see demos, type '.PBSfigs()'.
-----------------------------------------------------------


WARNING: RCall.jl Loading required package: maptools
Loading required package: sp
Checking rgeos availability: FALSE
 	Note: when rgeos is not available, polygon geometry 	computations in maptools depend on gpclib,
 	which has a restricted licence. It is disabled by default;
 	to enable gpclib, type gpclibPermit()
Loading required package: foreign

Now, perform the mapping.


In [5]:
R"ggplot($recorded, aes(map_id=FIPS)) +
geom_map(aes(fill=IR_To), map=shapes) +
expand_limits(x=c(-2500000, 2500000), y=c(-1.4e6, 1.6e6)) +
theme_bw() + theme(legend.justification=c(0,0), legend.position=c(0,0)) + xlab('') + ylab('')"


Out[5]:
RCall.RObject{RCall.VecSxp}

To compare this, we convert our results to USGS units, MGal/day.


In [6]:
simulated[:converted] = simulated[:totalirrigation] / (365 * 3785.41 / 1000)
recorded[:difference] = simulated[:converted] - recorded[:IR_To]
R"ggplot($recorded, aes(map_id=FIPS)) +
geom_map(aes(fill=difference), map=shapes) +
expand_limits(x=c(-2500000, 2500000), y=c(-1.4e6, 1.6e6)) +
theme_bw() + theme(legend.justification=c(0,0), legend.position=c(0,0)) + xlab('') + ylab('')"


Out[6]:
RCall.RObject{RCall.VecSxp}

Domestic Demand

Rather than a direct comparison of the population-driven domestic demand in the model and the USGS recorded values, we make a map of the apparent USGS use per person. Here is the domestic demand per person, in $1000 m^3 / year$.


In [7]:
m.external_parameters[:waterdemandperperson].value


Out[7]:
0.21001875

And below are the populations, in county FIPS order.


In [4]:
population =  m.external_parameters[:population].values[:, 1]


Out[4]:
3109-element Array{Float64,1}:
  54571.0
 182265.0
  27457.0
  22915.0
  57322.0
  10914.0
  20947.0
 118572.0
  34215.0
  25989.0
  43643.0
  13859.0
  25833.0
      ⋮  
  18106.0
  75450.0
   2484.0
  28205.0
   8667.0
  29116.0
  10247.0
  43806.0
  21294.0
  21118.0
   8533.0
   7208.0

We combine public supply and domestic self-supply to compare to the model's domestic demand.


In [7]:
recorded[:population] = population;
recorded[:percap] = (365 * 3785.41 / 1000) * (recorded[:DO_To] + recorded[:PS_To]) ./ recorded[:population]
R"ggplot($recorded, aes(map_id=FIPS)) +
geom_map(aes(fill=percap), map=shapes) + scale_fill_gradient(name='Per capita', trans='log', breaks=c(.0002, .002, .02, .2, 2, 20)) +
expand_limits(x=c(-2500000, 2500000), y=c(-1.4e6, 1.6e6)) +
theme_bw() + theme(legend.justification=c(0,0), legend.position=c(0,0)) + xlab('') + ylab('')"


Out[7]:
RCall.RObject{RCall.VecSxp}

Total water extraction

The total water in the model combines irrigation, domestic use, and exogenous thermodyanmic and livestock demands.


In [10]:
simulated = getdataframe(m, :WaterDemand, :totaldemand)


Out[10]:
regionstimetotaldemand
101001201519742.65928625
201003201538984.18938875
30100520159651.56833875
40100720156914.119496249999
501009201513393.6349475
60101120152485.7075175
70101320155643.59555625
801015201525372.424505
90101720157447.995821222293
100101920155845.30305375
110102120159414.71486625
120102320153076.56089625
130102520155618.97724875
140102720153354.584745
150102920153709.4066420778727
1601031201513446.011870081906
170103320151.7570777338049999e6
180103520153192.905625
190103720152492.5359562500003
2001039201511291.05665375
210104120155002.981742479571
2201043201519721.0812125
2301045201511203.470446250001
2401047201519684.101540880212
2501049201517830.247104415965
2601051201517545.415738294207
270105320158298.16310263306
28010552015181469.54894250003
290105720153897.45166875
300105920157695.37845
&vellip&vellip&vellip&vellip

In [11]:
R"ggplot($recorded, aes(map_id=FIPS)) +
geom_map(aes(fill=TO_To), map=shapes) +
expand_limits(x=c(-2500000, 2500000), y=c(-1.4e6, 1.6e6)) +
theme_bw() + theme(legend.justification=c(0,0), legend.position=c(0,0)) + xlab('') + ylab('')"


Out[11]:
RCall.RObject{RCall.VecSxp}

In [12]:
simulated[:converted] = simulated[:totaldemand] / (365 * 3785.41 / 1000)
recorded[:difference] = simulated[:converted] - recorded[:TO_To]
R"ggplot($recorded, aes(map_id=FIPS)) +
geom_map(aes(fill=difference), map=shapes) +
expand_limits(x=c(-2500000, 2500000), y=c(-1.4e6, 1.6e6)) +
theme_bw() + theme(legend.justification=c(0,0), legend.position=c(0,0)) + xlab('') + ylab('')"