compute_contributions_debug.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import rasterio as rs\n",
    "import geopandas as gpd\n",
    "from netCDF4 import Dataset\n",
    "import os\n",
    "import datetime\n",
    "import sys\n",
    "import gdal\n",
    "\n",
    "#reg = sys.argv[1] # pull the argument passed to the script\n",
    "reg = '01'\n",
    "fl = './data/nhrus/clean_AEA/nhru_%s_clean.shp'%reg # shapefile\n",
    "idxraster = './data/livneh_idx_laea.tiff' # the index raster to use with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cell(geom,gt=[],rb=[]):\n",
    "    ''' Grab index cell value for an hru that is too small to create a tiff for.\n",
    "    geom = hru geometry \n",
    "    gt = geotransform\n",
    "    rb = raster band to pick data from\n",
    "    '''\n",
    "    # grab the centrioid\n",
    "    mx = geom.centroid.x\n",
    "    my = geom.centroid.y\n",
    "    \n",
    "    # transform to array coordinates\n",
    "    px = int((mx-gt[0])/gt[1]) # x pixle\n",
    "    py = int((my-gt[3])/gt[5]) # y pixle\n",
    "    \n",
    "    # extract the value\n",
    "    intval = rb.ReadAsArray(px,py,1,1)\n",
    "    return intval[0][0]\n",
    "\n",
    "def process_tiffs(df,gt=[],rb=[]):\n",
    "    '''Process the tiffs or use the cell finder if the tiff does not exist\n",
    "    Inputs:\n",
    "    df = geopandas data frame of the cooresponding regional shapefile\n",
    "    gt = geotransform\n",
    "    rb = raster band\n",
    "    '''\n",
    "    \n",
    "    fl = './data/nhrus/AEA_tiffs/livneh_HUC_%s_nhruID_%s.tiff'%(df.region,df.hru_id_nat) # path to the tiffs\n",
    "\n",
    "    if os.path.isfile(fl) == True: # only proceed if the tiff exists\n",
    "        with rs.open(fl) as ds:\n",
    "            rast = ds.read(1)\n",
    "            nd = ds.NoData\n",
    "\n",
    "        n,m = rast.shape\n",
    "        rast.shape = n*m\n",
    "        rast = rast[rast!=nd] # remove no data cells\n",
    "        k = len(rast) # number of smaller cells in hru\n",
    "\n",
    "        cells = np.unique(rast)\n",
    "        #print(len(cells))\n",
    "        percents = []\n",
    "        for cell in cells:\n",
    "            percents.append(len(rast[rast==cell])/k) # divide by the total cells in the basin to get the propotion of each cell in the basin\n",
    "\n",
    "        cells = list(cells)\n",
    "        return cells,percents\n",
    "    \n",
    "    elif os.path.isfile(fl) == False: # if the raster does not exist, find the grid cell that the hru centroid occupies\n",
    "        print('W - National HRU %s tiff missing, finding index value at centroid'%df.hru_id_nat)\n",
    "        return [get_cell(df.geometry,gt = gt,rb = rb)],[1.]\n",
    "\n",
    "def compute_contributions(fl,idxraster = idxraster):\n",
    "    '''Compute grid cell contributions to hrus within a region'''\n",
    "    \n",
    "    reg = fl.split('_')[-2]\n",
    "    \n",
    "    print('S - Starting Region: %s!'%reg)\n",
    "    print('S - Using input shapefile: %s'%fl)\n",
    "    print('S - Using index raster: %s'%idxraster)\n",
    "    \n",
    "    # load the index raster and pull the geotransformation:\n",
    "    src_ds = gdal.Open(idxraster)\n",
    "    gt = src_ds.GetGeoTransform()\n",
    "    rb = src_ds.GetRasterBand(1)\n",
    "\n",
    "    dat = gpd.read_file(fl)\n",
    "    \n",
    "    cells,percents = zip(*dat.apply(process_tiffs,axis=1,gt=gt,rb=rb)) # run the aggregation function\n",
    "    \n",
    "    dat['reg'] = reg\n",
    "    \n",
    "    dat['cells'] = cells # insert results back into the dataframe\n",
    "    dat['percents'] = percents\n",
    "    \n",
    "    # remove geometry from the data frame:\n",
    "    del dat['geometry']\n",
    "    outputFile = './data/livneh_huc_%s_cell_contrib.pcl'%reg\n",
    "    print('S - Writing contribution data to: %s'%outputFile)\n",
    "    dat.to_pickle(outputFile)\n",
    "    print('S - %s Complete!'%reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S - Starting Region: 01!\n",
      "S - Using input shapefile: ./data/nhrus/clean_AEA/nhru_01_clean.shp\n",
      "S - Using index raster: ./data/livneh_idx_laea.tiff\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "(\"'rasterio._io.RasterReader' object has no attribute 'NoData'\", 'occurred at index 0')",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-e109267054fa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcompute_contributions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-2-0a82c7a8df1d>\u001b[0m in \u001b[0;36mcompute_contributions\u001b[0;34m(fl, idxraster)\u001b[0m\n\u001b[1;32m     66\u001b[0m     \u001b[0mdat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 68\u001b[0;31m     \u001b[0mcells\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpercents\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprocess_tiffs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mgt\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mrb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# run the aggregation function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     69\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     70\u001b[0m     \u001b[0mdat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'reg'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreg\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/py36/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36mapply\u001b[0;34m(self, func, axis, broadcast, raw, reduce, args, **kwds)\u001b[0m\n\u001b[1;32m   4875\u001b[0m                         \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4876\u001b[0m                         \u001b[0mreduce\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4877\u001b[0;31m                         ignore_failures=ignore_failures)\n\u001b[0m\u001b[1;32m   4878\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4879\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_apply_broadcast\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/py36/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_apply_standard\u001b[0;34m(self, func, axis, ignore_failures, reduce)\u001b[0m\n\u001b[1;32m   4971\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4972\u001b[0m                 \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries_gen\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4973\u001b[0;31m                     \u001b[0mresults\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4974\u001b[0m                     \u001b[0mkeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4975\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/py36/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36mf\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m   4852\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mkwds\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0margs\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mufunc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4853\u001b[0m             \u001b[0;32mdef\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4854\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4855\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4856\u001b[0m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-2-0a82c7a8df1d>\u001b[0m in \u001b[0;36mprocess_tiffs\u001b[0;34m(df, gt, rb)\u001b[0m\n\u001b[1;32m     30\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mrs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfl\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m             \u001b[0mrast\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m             \u001b[0mnd\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mNoData\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m         \u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrast\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: (\"'rasterio._io.RasterReader' object has no attribute 'NoData'\", 'occurred at index 0')"
     ]
    }
   ],
   "source": [
    "compute_contributions(fl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = gpd.read_file(fl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "geom = dat.loc[dat.hru_id_nat == 5].geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4    POINT (2188401.747535211 -69085.0446384551)\n",
       "dtype: object"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geom.centroid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rs.open(idxraster) as ds:\n",
    "    dat = ds.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(None,)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.nodata"
   ]
  },
  {
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}