THIS SCRIPT COMPARES WATERSHED AREAS ACCUMULATED IN STREAMCAT VS. AGAP

areaCompare.py

# -*- coding: utf-8 -*-
"""
Created on Mon Feb 06 16:02:50 2017

@author: Rdebbout
"""

import os
import h5py
import struct
import decimal
import datetime
import itertools
import numpy as np
import pandas as pd

def getUpstream(com, c, l ,u):
    itemindex = int(np.where(c == com)[0])
    n = l[:itemindex].sum()
    arrlen = l[itemindex]
    return u[n:n+arrlen]

def dbf2DF(f, upper=True):
    data = list(dbfreader(open(f, 'rb')))
    if upper == False:    
        return pd.DataFrame(data[2:], columns=data[0])
    else:
        return pd.DataFrame(data[2:], columns=map(str.upper,data[0]))    

def dbfreader(f):
    numrec, lenheader = struct.unpack('<xxxxLH22x', f.read(32))    
    numfields = (lenheader - 33) // 32
    fields = []
    for fieldno in xrange(numfields):
        name, typ, size, deci = struct.unpack('<11sc4xBB14x', f.read(32))
        name = name.replace('\0', '')       # eliminate NULs from string   
        fields.append((name, typ, size, deci))
    yield [field[0] for field in fields]
    yield [tuple(field[1:]) for field in fields]
    terminator = f.read(1)
    assert terminator == '\r'
    fields.insert(0, ('DeletionFlag', 'C', 1, 0))
    fmt = ''.join(['%ds' % fieldinfo[2] for fieldinfo in fields])
    fmtsiz = struct.calcsize(fmt)
    for i in xrange(numrec):
        record = struct.unpack(fmt, f.read(fmtsiz))
        if record[0] != ' ':
            continue                        # deleted record
        result = []
        for (name, typ, size, deci), value in itertools.izip(fields, record):
            if name == 'DeletionFlag':
                continue
            if typ == "N":
                value = value.replace('\0', '').lstrip()
                if value == '':
                    value = 0
                elif deci:
                    value = decimal.Decimal(value)
                else:
                    value = int(value)
            elif typ == 'C':
                value = value.rstrip()                                   
            elif typ == 'D':
                try:
                    y, m, d = int(value[:4]), int(value[4:6]), int(value[6:8])
                    value = datetime.date(y, m, d)
                except:
                    value = None
            elif typ == 'L':
                value = (value in 'YyTt' and 'T') or (value in 'NnFf' and 'F') or '?'
            elif typ == 'F':
                value = float(value)
            result.append(value)
        yield result
        
###############################################################################

if __name__ == '__main__':
    
    os.chdir('D:\\Projects\\CDI')
    inputs = np.load('./zoneInputs.npy').item()
    nhd = 'D:/NHDPlusV21'
    h = 'naqwaarea.h5'
    hdf = h5py.File(h,'r')
    sc = 'L:/Priv/CORFiles/Geospatial_Library/Data/Project/StreamCat/Allocation_and_Accumulation'
    cats = np.load('allCatCOMsZone.npy')[:,0]
    for zone in inputs.keys():
        print zone
        tbl = pd.read_csv('%s/precip_%s.csv' % (sc,zone))
        tbl = tbl[['COMID','WsAreaSqKm']]
        hr = inputs[zone]
        pre = "%s/NHDPlus%s/NHDPlus%s" % (nhd, hr, zone)
        fl = dbf2DF('%s/NHDSnapshot/Hydrography/NHDFlowline.dbf' % pre)
        flns = fl.ix[~fl.FLOWDIR.isin(['Uninitialized'])].COMID.values    
        coms = flns[np.in1d(flns,cats)]
        if zone == '04':
            f = 'L:/Priv/CORFiles/Geospatial_Library/Data/Project/StreamCat/FTP_Staging/StreamCat/Documentation/DataProcessingAndQualityAssurance/QA_Files/ProblemStreamsR04.csv'
            no = pd.read_csv(f)[['COMID']].values[:,0]
            np.in1d(coms,no)
            coms = coms[~np.in1d(coms,no)]
        print '{:,} COMIDs in {}'.format(len(coms),zone)
        d = {'COMID':[],'AreaM':[]}
        for x in coms:
            tmp = hdf.get(str(x))
            if tmp == None:
                print x
                continue
            d['COMID'].append(x)
            d['AreaM'].append(int(np.array(tmp.get('totArea'))))
        if not len(d['COMID']) == len(d['AreaM']):
            print 'Zonos Problemos: %s' % zone
        difs = pd.DataFrame(d)
        them = pd.merge(tbl,difs,on='COMID')                
        them.columns.tolist()
        them['HDFAreaSqKm'] = (them.AreaM * 1e-6).apply(lambda x: round(x,4))
        them['diff'] = abs(them.WsAreaSqKm - them.HDFAreaSqKm)
        them.sort_values('diff',ascending=False,inplace=True)
        these = them.ix[them['diff'] > 0.0011]
        them = pd.merge(them,fl.ix[fl.COMID.isin(them.COMID)][['COMID','FTYPE']],on='COMID')
        them.to_csv('./compareArea/compAll_%s.csv' % zone, index=False)