Here we compare the results of the simulation to recorded extractions from the USGS.
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")
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]:
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]:
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))";
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]:
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]:
In [7]:
m.external_parameters[:waterdemandperperson].value
Out[7]:
And below are the populations, in county FIPS order.
In [4]:
population = m.external_parameters[:population].values[:, 1]
Out[4]:
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]:
In [10]:
simulated = getdataframe(m, :WaterDemand, :totaldemand)
Out[10]:
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]:
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('')"
Out[12]: