Lab2_Assignment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 2 -- Assignment\n",
    "\n",
    "## In this Lab you will:\n",
    " - Get county average corn yield data for Indiana (or another state of interest)\n",
    " - Wrangle the data using Excel to delete columns of little or no interest (you could do it in Python but we need to save time)\n",
    " - Follow the example of Lecture 2 to estimate the corn yield trend over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring in the packages we have used before.\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Data\n",
    "To get some data for this exercise click on the link: <a href=\"https://quickstats.nass.usda.gov/\" target=\"_blank\">the USDA National Agricultural Statistics Service</a>. I'm using Indiana as an example below, but you are welcome to change it to Iowa, Illinois, Ohio, or name your state. There could, however, be some special cases that would have to be considered for different states.\n",
    "\n",
    "On the quick stats page you should select: Survey; Crops; Field Crops; Corn; Yield; Corn Grain - Yield, Measured in Bu/Acre; County; Indiana (leave unchecked the rest of the choices).\n",
    "\n",
    "Then click on the `get data` box at the bottom of the page. A spreadsheet should show up in the browser.\n",
    "\n",
    "Then click on the `Spreadsheet` towards the top right. This should download a `.csv` file somewhere on your computer or in the docker if you are running this on `mybinder.org`. It will have a nasty long name. Rename it to something convenient.\n",
    "\n",
    "If you are running on `mybinder.org` you will need to get this `.csv` onto your computer in order to use Excel.\n",
    "\n",
    "Use Excel to remove un-needed columns.\n",
    "\n",
    "Upload it to the directory `Lab2Data` in your fork of the \n",
    "<a href=\"https://github.com/jvkrogmeier/HEC-Lab-2\" target=\"_blank\">JVK github repo for HEC-Lab-2</a>.\n",
    "\n",
    "If you are already running in `mybinder.org` you will need to close and restart to have your `.csv` file available ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Name your file and put the name in here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replace Lab2Data/XXX.csv with Lab2Data/YourFileName.csv\n",
    "\n",
    "my_file_name = 'Lab2Data/XXX.csv'# FIX THIS!\n",
    "\n",
    "\n",
    "my_file_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just run this block which will open the wrangled file for reading \n",
    "# only and print the first few lines. What would you change in the \n",
    "# block below to see more or fewer lines?\n",
    "\n",
    "fin = open(my_file_name, \"r\")\n",
    "\n",
    "for i in range(20):\n",
    "    line = fin.readline()\n",
    "    print(line, end='')\n",
    "\n",
    "fin.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Remove some unwanted non-counties\n",
    "\n",
    "There are some names in the County field that correspond to combinations of counties or such. First we need to read the entire file to see what unique county names are present. There are probably some combined counties that we will want to elimate. Since we are concentrating only on one state we don't need to keep that information ...\n",
    "\n",
    "Run the code below (and suitably modify where requested) to remove the undesired non-counties.\n",
    "\n",
    "In the code block below you need to modify the integer indices corresponding to the `Year` column and the `County` column in your csv file. I've set the variables to something that won't work. You need to change it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the wrangled file.\n",
    "\n",
    "fin = open(my_file_name, \"r\")\n",
    "line = fin.readline() #Read first line to move past header \n",
    "\n",
    "# Initialize some empty sets.\n",
    "years = set()\n",
    "counties = set()\n",
    "\n",
    "yearindex = 59 # FIX THIS!\n",
    "countyindex = 137 # FIX THIS!\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for line in fin:\n",
    "    fields = line.split(',')\n",
    "    yr = fields[yearindex].replace('\"', '')\n",
    "    years.add(yr)\n",
    "    c = fields[countyindex].replace('\"', '')\n",
    "    counties.add(c)\n",
    "                      \n",
    "fin.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unique years, unique state_county_pairs\n",
    "The set variable `years` will contain the unique years for which we have data. Similarly the set variable `counties` contains the unique locations for which we have data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_unique_years = len(years)\n",
    "print(f'The number of unique years is {N_unique_years:2d}.')\n",
    "print(f'The set of years is {years}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_unique_counties = len(counties)\n",
    "print(f'The number of unique counties is {N_unique_counties:2d}.')\n",
    "print(f'The set of counties is {counties}.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bad County Names\n",
    "There are a certain number of county names in the `counties` list, which are not real county names. Define text strings that contain the bad county names so that we can ignore them later. By looking at the set of unique counties printed above, find the names we want to ignore and put them into the string variables below. For the Indiana data there are two bad names. Other states may have more or less ...\n",
    "\n",
    "`BCN1` and `BCN2` are string variables. Don't forget the quotes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bad county name 1, 2, 3 ... (BCN1, BCN2, etc)\n",
    "\n",
    "BCN1 = Blah # FIX THIS!\n",
    "BCN2 = BlahBlah # FIX THIS!\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just to check \n",
    "\n",
    "print(f'The number of unique counties (before removal) is {N_unique_counties:2d}.')\n",
    "\n",
    "# Remove the bad county names ...\n",
    "counties.remove(BCN1)\n",
    "counties.remove(BCN2)\n",
    "\n",
    "N_unique_counties = len(counties)\n",
    "\n",
    "print(f'The number of unique counties (after removal) is {N_unique_counties:2d}.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Convert the Sets to Lists and Sort\n",
    "We originally used sets to hold the `years` and `counties` data in order to impose uniqueness since sets never have duplicate elements. Therefore, when we were originally reading the data if the same year or same state-county pair was read multiple times (and they were), then we would automatically avoid adding duplicates.\n",
    "\n",
    "Finally, we convert the sets to lists in order to sort them. We particularly want the `years` variable ordered numerically increasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert them to lists. You need to replace XXX with the correct Python\n",
    "# command to convert the sets to lists.\n",
    "\n",
    "years = XXX(years) # FIX THIS!\n",
    "counties = XXX(counties) # FIX THIS!\n",
    "\n",
    "\n",
    "\n",
    "# Sort the years and counties lists. This is done so that we\n",
    "# can refer to the data in the numpy array and have the order there\n",
    "# correspond to the counties and year of interest.\n",
    "\n",
    "# Two commands need to be run here in order to sort the years \n",
    "# and counties lists ...\n",
    "\n",
    "put command 1 here # FIX THIS!\n",
    "put command 2 here # FIX THIS!\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Print just to see what the ordering looks like\n",
    "print(years)\n",
    "print()\n",
    "print(counties)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Create the `isdata` and `ydata` arrays\n",
    "Run the code blocks below (modify as needed) to create a numpy array that contains the information about whether a particular (county, year) has data in your file. The array will have rows indexed by county and columns indexed by year.\n",
    "\n",
    "You will also define an array `ydata` to contain the yield data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The variable \"isdata\" is initialized with zeros. When a one appears\n",
    "# at a particular row and column it indicates that the corresponding\n",
    "# county and year has yield data.\n",
    "\n",
    "isdata = np.zeros((len(counties), len(years)))\n",
    "\n",
    "# The variable \"ydata\" is also initialized with zeros.\n",
    "\n",
    "ydata = np.zeros((len(counties), len(years)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Populating data into `isdata` and `ydata`\n",
    "The next cell reads data from the `Lab2Data/XXX.csv` file and populates the numpy arrays `isdata` and `ydata`. Remember these two arrays were initialized as all zeros.\n",
    "\n",
    "Below we will use the `index()` method, which returns the numeric index of the first occurrence of an item in a list. We made our `counties` and `years` lists to have unique items."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the wrangled file.\n",
    "\n",
    "fin = open(my_file_name, \"r\")\n",
    "line = fin.readline() #Read first line to move past header \n",
    "\n",
    "# The year values are located in the column indexed by yearindex\n",
    "# which you set in a previous cell. = County values in the fourth \n",
    "# are indexed by countyindex, also set previously. You now need to\n",
    "# set the variable yieldindex to refer to the column of the csv\n",
    "# file which contains County average yield values.\n",
    "\n",
    "yieldindex = 11 # FIX THIS!\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Note how we use the python method .index in order to properly\n",
    "# align the data values in the numpy arrays with the list values\n",
    "# in years and state_county_pairs.\n",
    "\n",
    "for line in fin:\n",
    "    fields = line.split(',') #Separate line into individual items\n",
    "    yr = fields[yearindex].replace('\"', '') #Delete extra double quotes\n",
    "    c = fields[countyindex].replace('\"', '') #Delete extra double quotes\n",
    "    if c in counties:\n",
    "        rindex = counties.index(c)\n",
    "        cindex = years.index(yr)\n",
    "        isdata[rindex, cindex] = 1 #Insert \"1\" into isdata in proper place\n",
    "        yd = fields[yieldindex]\n",
    "        yd = yd.replace('\\n', '') #Delete newline character at the end\n",
    "        yd = yd.replace('\"', '') #Delete extra double quotes\n",
    "        ydata[rindex, cindex] = float(yd) #Insert a numeric yield value in the proper place\n",
    "                      \n",
    "fin.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. The `if c in counties:` block was not used in the Lecture 3 code. Why is it needed here and not there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just to look at the array\n",
    "print(ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just to look at the array\n",
    "print(isdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just run this code, which prepares a meshgrid for contour plots ...\n",
    "iyrs = range(len(isdata[0,:]))\n",
    "iscprs = range(len(isdata[:,0]))\n",
    "IYRS, ISCPRS = np.meshgrid(iyrs,iscprs)\n",
    "\n",
    "# Defining a function to evaluate values from isdata\n",
    "def ID(x,y):\n",
    "    return isdata[x,y]\n",
    "\n",
    "# Evaluate on the mesh\n",
    "Z = ID(ISCPRS,IYRS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make filled contour plot and label axes, etc. ...\n",
    "plt.contourf(IYRS,ISCPRS,Z, cmap='winter');\n",
    "plt.ylabel('Indices of counties')\n",
    "plt.xlabel('Indices of years (starting from 1929)')\n",
    "plt.title('Data Available (green = have data, blue = no data)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat only for the actual county average yield values\n",
    "def Y(x,y):\n",
    "    return ydata[x,y]\n",
    "Z = Y(ISCPRS,IYRS)\n",
    "\n",
    "plt.contourf(IYRS,ISCPRS,Z, cmap='YlOrRd');\n",
    "plt.ylabel('Indices of counties')\n",
    "plt.xlabel('Indices of years (starting from 1929)')\n",
    "plt.title('County average yield (bushels/acre)')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the fraction of complete data in this data set ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frac = (len(isdata[0,:])*len(isdata[:,0]) - sum(sum(isdata)))/(len(isdata[0,:])*len(isdata[:,0]))\n",
    "\n",
    "print(f'Fraction of unavailable data is {frac:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. What is the first thing you notice when comparing the situation with Indiana corn yield data with wheat yield data as in the Lect 2 code?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Some of the Indiana counties have complete yield data (i.e., for every one of the 92 years). Run the code blocks below to find out which counties they are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize a vector to hold info about which counties contain\n",
    "# data for every year. Then sequentially multiply all\n",
    "# of the columns together.\n",
    "all = np.ones(len(counties))\n",
    "\n",
    "for i in range(len(years)):\n",
    "    all = all*isdata[:,i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the result\n",
    "print(all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clearly some state-county pairs have data for all years. How many?\n",
    "all.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result of multiplying all rows of `isdata` ...\n",
    "According to this test there are 45 Indiana counties which have data in all 92 years. Which counties are they?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c92 = []\n",
    "i = 0\n",
    "for c in counties:\n",
    "    if all[i] > 0:\n",
    "        c92.append(c)\n",
    "    i = i + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(c92)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's use these `c92` counties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting the data for the restricted state_county pairs set ...\n",
    "In the code block below we want to create a list containing the indices from the entire collection of Indiana counties which correspond to the `c92`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rindex = []\n",
    "for c in counties:\n",
    "    if c in c92:\n",
    "        rindex.append(counties.index(c))\n",
    "\n",
    "# This creates the column index set corresponding to the years \n",
    "# for which data is available (which is actually every year)\n",
    "cindex = []\n",
    "for y in years:\n",
    "    cindex.append(years.index(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Select a smaller array with continguous data ...\n",
    "Out of the large array of size 92 x 92 containing county average yield data over various years we want to extract the smaller 45 x 92 array containing county average yield data for our selected state, county pairs of interest and our selected years of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the 45 rows corresponding to the counties of interest\n",
    "\n",
    "select_ydata = XXX[YYY,ZZZ] # FIX THIS!\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the individual county average yields on the same axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.style.use('classic')\n",
    "\n",
    "for k in range(select_ydata.shape[0]):\n",
    "    plt.plot(select_ydata[k,:])\n",
    "    \n",
    "plt.title(\"County Average Yields vs. Year\")\n",
    "plt.xlabel(\"Year since 1929\")\n",
    "plt.ylabel(\"Yield (bu/acre)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. What happened in the year just past 80 in the plot above?\n",
    "\n",
    "## 11. Compute the mean over all the counties in the array `select_ydata`. \n",
    "\n",
    "The resulting mean will still depend on the year, which will allow us to see trend in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "select_ydata_mean = np.XXX(select_ydata,axis=YYY) # FIX THIS!\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.style.use('classic')\n",
    "plt.plot(select_ydata_mean,'*')\n",
    "plt.title(\"Mean Yields vs. Year\")\n",
    "plt.xlabel(\"Year since 1929\")\n",
    "plt.ylabel(\"Yield (bu/acre)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Straight Line Fitting via Least Squares ...\n",
    "\n",
    "<img align=\"left\" src='Data2/LS-notes-p1.png' width=\"300\"/>\n",
    "<img align=\"left\" src='Data2/LS-notes-p2.png' width=\"300\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First lets fit a line to the mean over the counties of \n",
    "# interest and see how that varies over time. The items below \n",
    "# correspond directly to the entries in the LS notes above.\n",
    "\n",
    "y = select_ydata_mean\n",
    "N = len(y)\n",
    "x = range(N) # Remember we are using years starting from 0 in 1929\n",
    "xy = np.multiply(x,y)\n",
    "xx = np.multiply(x,x)\n",
    "xsum = sum(x)\n",
    "ysum = sum(y)\n",
    "xxsum = sum(xx)\n",
    "xysum = sum(xy)\n",
    "\n",
    "A = np.zeros((2,2))\n",
    "A[0,0] = N\n",
    "A[0,1] = xsum\n",
    "A[1,0] = xsum\n",
    "A[1,1] = xxsum\n",
    "\n",
    "b = np.zeros(2)\n",
    "b[0] = ysum\n",
    "b[1] = xysum\n",
    "\n",
    "z = np.linalg.solve(A, b)\n",
    "print(z)\n",
    "\n",
    "# Now plot the best fit line on top of the data\n",
    "fig = plt.figure()\n",
    "plt.style.use('classic')\n",
    "plt.plot(x,y,'*')\n",
    "plt.title(\"Mean of Yields vs. Year\")\n",
    "plt.xlabel(\"Year since 1929\")\n",
    "plt.ylabel(\"Yield (bu/acre)\")\n",
    "plt.grid()\n",
    "yy = z[0] + z[1]*x\n",
    "plt.plot(x,yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. What is the Yield Improvement Trend in bu/acre per year?\n",
    "\n",
    "There is quite a difference between corn and wheat breeding!\n",
    "\n",
    "## 13. Google to see what you can find on the corn yield improvement trend. Does this data agree with what you find?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# With a for loop ...\n",
    "fig = plt.figure()\n",
    "plt.style.use('classic')\n",
    "plt.title(\"Mean of Yields vs. Year\")\n",
    "plt.xlabel(\"Year since 1970\")\n",
    "plt.ylabel(\"Yield (bu/acre)\")\n",
    "plt.grid()\n",
    "\n",
    "for k in range(28):\n",
    "    y = select_ydata[k,:]\n",
    "    N = len(y)\n",
    "    x = range(N) # Remember we are years starting from 0\n",
    "    xy = np.multiply(x,y)\n",
    "    xx = np.multiply(x,x)\n",
    "    xsum = sum(x)\n",
    "    ysum = sum(y)\n",
    "    xxsum = sum(xx)\n",
    "    xysum = sum(xy)\n",
    "\n",
    "    A = np.zeros((2,2))\n",
    "    A[0,0] = N\n",
    "    A[0,1] = xsum\n",
    "    A[1,0] = xsum\n",
    "    A[1,1] = xxsum\n",
    "\n",
    "    b = np.zeros(2)\n",
    "    b[0] = ysum\n",
    "    b[1] = xysum\n",
    "\n",
    "    z = np.linalg.solve(A, b)\n",
    "\n",
    "    # Now plot the best fit line on top of the data\n",
    "    plt.plot(x,y,'.')\n",
    "    yy = z[0] + z[1]*x\n",
    "    plt.plot(x,yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}