#rootname = "/media/p/5F5B-8FCB/PROJECTS/UMAR/Data/chem/" #thumb on ubuntu
rootname = "E:\\PROJECTS\\UMAR\\Data\\chem\\" #thumb on windows

WQPResultsFile = rootname + "result.csv"
WQPStationFile = rootname + "station.csv"
SDWISFile = rootname + "SDWIS_Cache.txt"
AGStationsFile = rootname + "AG_Stations_Cache.csv"
AGResultsFile = rootname + "AG_byparam.csv"
UGSFile = rootname + "UGS_Cache.txt"
STORLegStatFile = rootname + "UT_Cache_sta_001.txt"
STORLegResFile = rootname + "UT_Cache_res_001.txt"
STORParamFile = rootname + "parameter.txt"

fielddata = rootname + "FieldData.xlsx"
statelabresults0 = rootname + "utgs1402.txt"
statelabresults1 = rootname + "utgs1403.txt"
statelabresults2 = rootname + "utgs1501.txt"
statelabstations = rootname + "UtahStateLabStations.xlsx"

df = wc.WQP.WQPimportRes(WQPResultsFile)
df = wc.WQP.WQPmassageResults(df)
df

"""
Explanation: Import and Standardize Data
End of explanation
"""


Rdtypes = {"OrganizationIdentifier":np.str_, "OrganizationFormalName":np.str_, "ActivityIdentifier":np.str_, 
           "ActivityStartTime/Time":np.str_,
           "ActivityTypeCode":np.str_, "ActivityMediaName":np.str_, "ActivityMediaSubdivisionName":np.str_, 
           "ActivityStartDate":np.str_, "ActivityStartTime/Time":np.str_, "ActivityStartTime/TimeZoneCode":np.str_, 
           "ActivityEndDate":np.str_, "ActivityEndTime/Time":np.str_, "ActivityEndTime/TimeZoneCode":np.str_, 
           "ActivityDepthHeightMeasure/MeasureValue":np.float16, "ActivityDepthHeightMeasure/MeasureUnitCode":np.str_, 
           "ActivityDepthAltitudeReferencePointText":np.str_, "ActivityTopDepthHeightMeasure/MeasureValue":np.float16, 
           "ActivityTopDepthHeightMeasure/MeasureUnitCode":np.str_, 
           "ActivityBottomDepthHeightMeasure/MeasureValue":np.float16, 
           "ActivityBottomDepthHeightMeasure/MeasureUnitCode":np.str_, 
           "ProjectIdentifier":np.str_, "ActivityConductingOrganizationText":np.str_, 
           "MonitoringLocationIdentifier":np.str_, "ActivityCommentText":np.str_, 
           "SampleAquifer":np.str_, "HydrologicCondition":np.str_, "HydrologicEvent":np.str_, 
           "SampleCollectionMethod/MethodIdentifier":np.str_, "SampleCollectionMethod/MethodIdentifierContext":np.str_, 
           "SampleCollectionMethod/MethodName":np.str_, "SampleCollectionEquipmentName":np.str_, 
           "ResultDetectionConditionText":np.str_, "CharacteristicName":np.str_, "ResultSampleFractionText":np.str_, 
           "ResultMeasureValue":np.str_, "ResultMeasure/MeasureUnitCode":np.str_, "MeasureQualifierCode":np.str_, 
           "ResultStatusIdentifier":np.str_, "StatisticalBaseCode":np.str_, "ResultValueTypeName":np.str_, 
           "ResultWeightBasisText":np.str_, "ResultTimeBasisText":np.str_, "ResultTemperatureBasisText":np.str_, 
           "ResultParticleSizeBasisText":np.str_, "PrecisionValue":np.str_, "ResultCommentText":np.str_, 
           "USGSPCode":np.str_, "ResultDepthHeightMeasure/MeasureValue":np.float16, 
           "ResultDepthHeightMeasure/MeasureUnitCode":np.str_, "ResultDepthAltitudeReferencePointText":np.str_, 
           "SubjectTaxonomicName":np.str_, "SampleTissueAnatomyName":np.str_, 
           "ResultAnalyticalMethod/MethodIdentifier":np.str_, "ResultAnalyticalMethod/MethodIdentifierContext":np.str_, 
           "ResultAnalyticalMethod/MethodName":np.str_, "MethodDescriptionText":np.str_, "LaboratoryName":np.str_, 
           "AnalysisStartDate":np.str_, "ResultLaboratoryCommentText":np.str_, 
           "DetectionQuantitationLimitTypeName":np.str_, "DetectionQuantitationLimitMeasure/MeasureValue":np.str_, 
           "DetectionQuantitationLimitMeasure/MeasureUnitCode":np.str_, "PreparationStartDate":np.str_, 
           "ProviderName":np.str_} 

dt = [6,56,61]

"""
Explanation: WQP
The following html addresses are REST-based queries to download WQP data from the <a href="http://www.waterqualitydata.us/portal/">WQP portal</a>. If you click on them, they will produce zipped csv files that can be opened and processed with the code below. Originally, the code directly applied these links, but the files are large and take a lot of time to download.
Station data address:
http://waterqualitydata.us/Station/search?statecode=US%3A49&countycode=US%3A49%3A005&sampleMedia=Water&characteristicType=Information%3BInorganics%2C+Major%2C+Metals%3BInorganics%2C+Major%2C+Non-metals%3BInorganics%2C+Minor%2C+Metals%3BInorganics%2C+Minor%2C+Non-metals%3BNot+Assigned%3BNutrient%3BPhysical%3BStable+Isotopes&mimeType=csv&zip=yes&sorted=no
Result data address:
http://waterqualitydata.us/Result/search?statecode=US%3A49&countycode=US%3A49%3A005&sampleMedia=Water&characteristicType=Information%3BInorganics%2C+Major%2C+Metals%3BInorganics%2C+Major%2C+Non-metals%3BInorganics%2C+Minor%2C+Metals%3BInorganics%2C+Minor%2C+Non-metals%3BNot+Assigned%3BNutrient%3BPhysical%3BStable+Isotopes&mimeType=csv&zip=yes&sorted=no
WQP Results
Define data type of each field in the WQP database.  This allows for easy import of data. Everything under this header can be acheived using wc.WQP.WQPimportRes(WQPResultsFile) then wc.WQP.WQPmassageResults(df)
End of explanation
"""


WQP = pd.read_csv(WQPResultsFile, dtype=Rdtypes, parse_dates=dt)

"""
Explanation: Read csv data into python.
End of explanation
"""


ResFieldDict = {"AnalysisStartDate":"AnalysisDate", "ResultAnalyticalMethod/MethodIdentifier":"AnalytMeth", 
                "ResultAnalyticalMethod/MethodName":"AnalytMethId", "ResultDetectionConditionText":"DetectCond", 
                "ResultLaboratoryCommentText":"LabComments", "LaboratoryName":"LabName", 
                "DetectionQuantitationLimitTypeName":"LimitType", "DetectionQuantitationLimitMeasure/MeasureValue":"MDL", 
                "DetectionQuantitationLimitMeasure/MeasureUnitCode":"MDLUnit", "MethodDescriptionText":"MethodDescript", 
                "OrganizationIdentifier":"OrgId", "OrganizationFormalName":"OrgName", "CharacteristicName":"Param", 
                "ProjectIdentifier":"ProjectId", "MeasureQualifierCode":"QualCode", "ResultCommentText":"ResultComment", 
                "ResultStatusIdentifier":"ResultStatus", "ResultMeasureValue":"ResultValue", 
                "ActivityCommentText":"SampComment", "ActivityDepthHeightMeasure/MeasureValue":"SampDepth", 
                "ActivityDepthAltitudeReferencePointText":"SampDepthRef", 
                "ActivityDepthHeightMeasure/MeasureUnitCode":"SampDepthU", "SampleCollectionEquipmentName":"SampEquip", 
                "ResultSampleFractionText":"SampFrac", "ActivityStartDate":"SampleDate", "ActivityIdentifier":"SampleId", 
                "ActivityStartTime/Time":"SampleTime", "ActivityMediaSubdivisionName":"SampMedia", 
                "SampleCollectionMethod/MethodIdentifier":"SampMeth", "SampleCollectionMethod/MethodName":"SampMethName", 
                "ActivityTypeCode":"SampType", "MonitoringLocationIdentifier":"StationId", 
                "ResultMeasure/MeasureUnitCode":"Unit", "USGSPCode":"USGSPCode",
                "ActivityStartDate":"StartDate","ActivityStartTime/Time":"StartTime"}

WQP.rename(columns=ResFieldDict,inplace=True)

"""
Explanation: Rename columns to match with other data later.
End of explanation
"""


resdroplist = ["ActivityBottomDepthHeightMeasure/MeasureUnitCode", "ActivityBottomDepthHeightMeasure/MeasureValue", 
               "ActivityConductingOrganizationText", "ActivityEndDate", "ActivityEndTime/Time", 
               "ActivityEndTime/TimeZoneCode", "ActivityMediaName", "ActivityStartTime/TimeZoneCode", 
               "ActivityTopDepthHeightMeasure/MeasureUnitCode", "ActivityTopDepthHeightMeasure/MeasureValue", 
               "HydrologicCondition", "HydrologicEvent", "PrecisionValue", "PreparationStartDate", "ProviderName", 
               "ResultAnalyticalMethod/MethodIdentifierContext", "ResultDepthAltitudeReferencePointText", 
               "ResultDepthHeightMeasure/MeasureUnitCode", "ResultDepthHeightMeasure/MeasureValue", 
               "ResultParticleSizeBasisText", "ResultTemperatureBasisText", 
               "ResultTimeBasisText", "ResultValueTypeName", "ResultWeightBasisText", "SampleAquifer", 
               "SampleCollectionMethod/MethodIdentifierContext", "SampleTissueAnatomyName", "StatisticalBaseCode", 
               "SubjectTaxonomicName","StartTime","StartDate","StartTime","StartDate"]

"""
Explanation: Define unneeded columns that will be dropped to save memory.
End of explanation
"""


def datetimefix(x,format):
    '''
    This script cleans date-time errors
    
    input
    x = date-time string
    format = format of date-time string
    
    output 
    formatted datetime type
    '''
    d = str(x[0]).lstrip().rstrip()[0:10]
    t = str(x[1]).lstrip().rstrip()[0:5].zfill(5)
    try:
        int(d[0:2])
    except(ValueError,TypeError,NameError):
        return np.nan
    try:
        int(t[0:2])
        int(t[3:5])
    except(ValueError,TypeError,NameError):
        t = "00:00"
   
    if int(t[0:2])>23:
        t = "00:00"
    elif int(t[3:5])>59:
        t = "00:00"
    else:
        t = t[0:2].zfill(2) + ":" + t[3:5]
    return datetime.datetime.strptime(d + " " + t, format)

WQP["SampleDate"] = WQP[["StartDate","StartTime"]].apply(lambda x: datetimefix(x,"%Y-%m-%d %H:%M"),1)

"""
Explanation: Define a function to fix funky dates found in the WQP database.  This includes impossible dates or dates with too many numbers.
End of explanation
"""


WQP.drop(resdroplist,inplace=True,axis=1)

"""
Explanation: Drop unwanted fields
End of explanation
"""


WQP['ResultValue'] = WQP['ResultValue'].convert_objects(convert_numeric=True)
WQP['MDL'] = WQP['MDL'].convert_objects(convert_numeric=True)

"""
Explanation: Convert result values and the MDL values to numeric fields from string fields.
End of explanation
"""


WQP['StationId'] = WQP['StationId'].str.replace('_WQX-','-')

"""
Explanation: Remove station duplicates by removing the preceding 'WQX' found in the station id field.
End of explanation
"""


#standardize all ug/l data to mg/l
def unitfix(x):
    z = str(x).lower()
    if z == "ug/l":
        return "mg/l"
    elif z == "mg/l":
        return "mg/l"
    else:
        return x

WQP.Unit = WQP.Unit.apply(lambda x: str(x).rstrip(), 1)
WQP.ResultValue = WQP[["ResultValue","Unit"]].apply(lambda x: x[0]/1000 if str(x[1]).lower()=="ug/l" else x[0], 1)
WQP.Unit = WQP.Unit.apply(lambda x: unitfix(x),1)

"""
Explanation: Standardize all ug/l data to mg/l by changing out the unit descriptor and dividing by 1000.
End of explanation
"""


def parnorm(x):
    p = str(x[0]).rstrip().lstrip().lower()
    u = str(x[2]).rstrip().lstrip().lower()
    if p == 'nitrate' and u == 'mg/l as n':
        return 'Nitrate', x[1]*4.427, 'mg/l'
    elif p == 'nitrite' and u == 'mg/l as n':
        return 'Nitrite', x[1]*3.285, 'mg/l'
    elif p == 'ammonia-nitrogen' or p == 'ammonia-nitrogen as n' or p == 'ammonia and ammonium':
        return 'Ammonium', x[1]*1.288, 'mg/l'
    elif p == 'ammonium' and u == 'mg/l as n':
        return 'Ammonium', x[1]*1.288, 'mg/l'
    elif p == 'sulfate as s':
        return 'Sulfate', x[1]*2.996, 'mg/l'
    elif p in ('phosphate-phosphorus', 'phosphate-phosphorus as p','orthophosphate as p'):
        return 'Phosphate', x[1]*3.066, 'mg/l'
    elif (p == 'phosphate' or p == 'orthophosphate') and u == 'mg/l as p':
        return 'Phosphate', x[1]*3.066, 'mg/l'
    elif u == 'ug/l':
        return x[0], x[1]/1000, 'mg/l'
    else:
        return x[0], x[1], str(x[2]).rstrip()

WQP['Param'], WQP['ResultValue'], WQP['Unit'] = zip(*WQP[['Param','ResultValue','Unit']].apply(lambda x: parnorm(x),1))

"""
Explanation: Normalize common nutrients so that they are all in the same type of units.  For example, sometimes nitrate is reported "as nitrogen" and sometimes it is reported "as nitrate".  The difference between the two types of reporting is a factor of 4.427!
End of explanation
"""


WQPStat = pd.read_csv(WQPStationFile)

"""
Explanation: WQP Stations
Read in WQP station data.
End of explanation
"""


StatFieldDict = {"MonitoringLocationIdentifier":"StationId", "AquiferName":"Aquifer", "AquiferTypeName":"AquiferType", 
             "ConstructionDateText":"ConstDate", "CountyCode":"CountyCode", "WellDepthMeasure/MeasureValue":"Depth", 
             "WellDepthMeasure/MeasureUnitCode":"DepthUnit", "VerticalMeasure/MeasureValue":"Elev", 
             "VerticalAccuracyMeasure/MeasureValue":"ElevAcc", "VerticalAccuracyMeasure/MeasureUnitCode":"ElevAccUnit", 
             "VerticalCollectionMethodName":"ElevMeth", "VerticalCoordinateReferenceSystemDatumName":"ElevRef", 
             "VerticalMeasure/MeasureUnitCode":"ElevUnit", "FormationTypeText":"FmType", 
             "WellHoleDepthMeasure/MeasureValue":"HoleDepth", "WellHoleDepthMeasure/MeasureUnitCode":"HoleDUnit", 
             "HorizontalAccuracyMeasure/MeasureValue":"HorAcc", "HorizontalAccuracyMeasure/MeasureUnitCode":"HorAccUnit", 
             "HorizontalCollectionMethodName":"HorCollMeth", "HorizontalCoordinateReferenceSystemDatumName":"HorRef", 
             "HUCEightDigitCode":"HUC8", "LatitudeMeasure":"Lat_Y", "LongitudeMeasure":"Lon_X", 
             "OrganizationIdentifier":"OrgId", "OrganizationFormalName":"OrgName", "StateCode":"StateCode", 
             "MonitoringLocationDescriptionText":"StationComment", "MonitoringLocationName":"StationName", 
             "MonitoringLocationTypeName":"StationType"}

WQPStat.rename(columns=StatFieldDict,inplace=True)

"""
Explanation: Rename and simplify station data column names for later compilation.
End of explanation
"""


statdroplist = ["ContributingDrainageAreaMeasure/MeasureUnitCode", "ContributingDrainageAreaMeasure/MeasureValue", 
                "DrainageAreaMeasure/MeasureUnitCode", "DrainageAreaMeasure/MeasureValue", "CountryCode", "ProviderName", 
                "SourceMapScaleNumeric"]
WQPStat.drop(statdroplist,inplace=True,axis=1)

"""
Explanation: Define the fields to drop to save memory.
End of explanation
"""


TypeDict = {"Stream: Canal":"Stream", "River/Stream":"Stream", 
            "Stream: Canal":"Stream", "Well: Test hole not completed as a well":"Well"}
WQPStat.StationType = WQPStat["StationType"].apply(lambda x: TypeDict.get(x,x),1)
WQPStat.Elev = WQPStat.Elev.apply(lambda x: np.nan if x==0.0 else round(x,1), 1)

"""
Explanation: Make station types in the StationType field consistent for easier summary and compilation later on.
End of explanation
"""


WQPStat['StationId'] = WQPStat['StationId'].str.replace('_WQX-','-')
WQPStat.drop_duplicates(subset=['StationId'],inplace=True)

"""
Explanation: Remove preceding WQX from StationId field to remove duplicate station data created by legacy database.
End of explanation
"""


SDWIS = pd.read_csv(SDWISFile)

def sampid(x):
    return "SDWIS" + str(x[0]) + str(x[1]) + str(x[2])[:-7]

def statid(x):
    return "SDWIS" + str(x[0]) + str(x[1]) 

def statnm(x):
    return str(str(x[0]) + " " + str(x[1])).title()

SDWIS["StationId"] = SDWIS[["FED_NM","ST_ID"]].apply(lambda x: statid(x),1)
SDWIS["StationName"] = SDWIS[["SYS_NM","FAC_NM"]].apply(lambda x: statnm(x),1)
SDWIS["SampleId"] = SDWIS[["FED_NM","ST_ID","DTE"]].apply(lambda x: sampid(x),1)
SDWIS["OrgId"] = "UDDW"
SDWIS["OrgName"] = "Utah Division of Drinking Water"
SDWIS["Elev"] = SDWIS["Z"].apply(lambda x: round(x*3.2808,1),1)
SDWIS["Unit"] = SDWIS["UOM_CODE"].apply(lambda x: str(x).lower(),1)
SDWIS["MDLUnit"] = SDWIS["DET_UOM"].apply(lambda x: str(x).lower(),1)
SDWIS["Param"] = SDWIS["ANLY_NM"].apply(lambda x: str(x).title().rstrip(),1)

SDWISFields ={"DTE":"SampleDate", "TYPE_CODE":"StationType",
              "CONC_MSR":"ResultValue", "DET_LIM":"MDL",
              "Y":"Lat_Y", "X":"Lon_X"}  
SDWIS.rename(columns=SDWISFields,inplace=True)

def datetimefixSDWIS(x,format):
    d = str(x).lstrip().rstrip()
    try:
        return datetime.datetime.strptime(d, "%m/%d/%Y %H:%M:%S")
    except(ValueError):
        return datetime.datetime.strptime(d, "%Y-%m-%d %H:%M:%S")

SDWIS["SampleDate"] = SDWIS["SampleDate"].apply(lambda x: datetimefixSDWIS(x,"%m/%d/%Y %H:%M:%S"),1)

"""
Explanation: SDWIS
SDWIS data were extracted from the Utah SDWIS database into ArcGIS 10.3.2 using the following SQL query.  NED 10m elevation and UTM coordinates were appended using ArcGIS.
SQL
SELECT UTV80.TINWSF.EXTERNAL_SYS_NUM AS "FED_NM",  UTV80.TINWSF.ST_ASGN_IDENT_CD AS "ST_ID",  UTV80.TINWSF.TYPE_CODE,  UTV80.TINWSYS.NAME AS "SYS_NM", UTV80.TINWSYS.D_PRIN_CNTY_SVD_NM AS "COUNTY", UTV80.TINWSF.NAME AS "FAC_NM",  UTV80.TINWSF.TINWSYS_IS_NUMBER AS "SY_NBR", UTV80.TINLOC.LATITUDE_MEASURE AS "Y",  UTV80.TINLOC.LONGITUDE_MEASURE AS "X",  UTV80.TINLOC.VERTICAL_MEASURE AS "Z", UTV80.TSASAMPL.COLLLECTION_END_DT AS "DTE", UTV80.TSAANLYT.NAME AS "ANLY_NM",  UTV80.TSASAR.CONCENTRATION_MSR AS "CONC_MSR", UTV80.TSASAR.TSASAR_IS_NUMBER AS "ID_NUM",  UTV80.TSASAR.UOM_CODE,   UTV80.TSASAR.DETECTN_LIMIT_NUM AS "DET_LIM", UTV80.TSASAR.DETECTN_LIM_UOM_CD AS "DET_UOM"  FROM  UTV80.TINWSF  INNER JOIN  UTV80.TINWSYS ON UTV80.TINWSF.TINWSYS_IS_NUMBER = UTV80.TINWSYS.TINWSYS_IS_NUMBER  INNER JOIN  UTV80.TINLOC ON UTV80.TINWSF.TINWSF_IS_NUMBER = UTV80.TINLOC.TINWSF_IS_NUMBER  INNER JOIN  UTV80.TSASMPPT ON UTV80.TINWSF.TINWSF_IS_NUMBER = UTV80.TSASMPPT.TINWSF0IS_NUMBER  INNER JOIN  UTV80.TSASAMPL ON UTV80.TSASMPPT.TSASMPPT_IS_NUMBER = UTV80.TSASAMPL.TSASMPPT_IS_NUMBER  INNER JOIN  UTV80.TSASAR ON UTV80.TSASAMPL.TSASAMPL_IS_NUMBER = UTV80.TSASAR.TSASAMPL_IS_NUMBER  INNER JOIN  UTV80.TSAANLYT ON UTV80.TSASAR.TSAANLYT_IS_NUMBER = UTV80.TSAANLYT.TSAANLYT_IS_NUMBER WHERE  (UTV80.TINWSYS.D_PRIN_CNTY_SVD_NM LIKE '%CACHE COUNTY%')   AND (UTV80.TSAANLYT.NAME LIKE '%NITRATE%' OR UTV80.TSAANLYT.NAME LIKE '%NITRITE%' OR UTV80.TSAANLYT.NAME LIKE '%AMMONI%' OR UTV80.TSAANLYT.NAME LIKE '%SULFATE%' OR UTV80.TSAANLYT.NAME LIKE '%TDS%' OR UTV80.TSAANLYT.NAME LIKE '%SODIUM%' OR UTV80.TSAANLYT.NAME LIKE '%FLUORIDE%' OR UTV80.TSAANLYT.NAME LIKE '%MAGNESIUM%' OR UTV80.TSAANLYT.NAME LIKE '%SELENIUM%' OR  UTV80.TSAANLYT.NAME LIKE '%CALCIUM%' OR UTV80.TSAANLYT.NAME LIKE '%CHLORIDE%' OR UTV80.TSAANLYT.NAME LIKE '%POTASSIUM%' OR UTV80.TSAANLYT.NAME LIKE '%SELENIUM%' OR UTV80.TSAANLYT.NAME LIKE '%SILICA%' OR UTV80.TSAANLYT.NAME LIKE '%IRON %' OR UTV80.TSAANLYT.NAME LIKE '%ALKA %' OR UTV80.TSAANLYT.NAME LIKE '%CONDUCTIVITY%' OR UTV80.TSAANLYT.NAME LIKE '%PH %' OR UTV80.TSAANLYT.NAME LIKE '%TEMP%' OR UTV80.TSAANLYT.NAME LIKE '%ARSENIC%' OR UTV80.TSAANLYT.NAME LIKE '%CARBON%' OR UTV80.TSAANLYT.NAME LIKE '%TRITIUM%' OR UTV80.TSAANLYT.NAME LIKE '%COPPER%' OR UTV80.TSAANLYT.NAME LIKE '%LEAD%' OR UTV80.TSAANLYT.NAME LIKE '%NITROGEN%' OR UTV80.TSAANLYT.NAME LIKE '%PHOSPHATE%' OR UTV80.TSAANLYT.NAME LIKE '%TDS%' OR UTV80.TSAANLYT.NAME LIKE '%ZINC%' OR UTV80.TSAANLYT.NAME LIKE '%IRON%' OR UTV80.TSAANLYT.NAME LIKE '%CHROMIUM%' ) ORDER BY UTV80.TINWSF.ST_ASGN_IDENT_CD
SQL
SELECT UTV80.TINWSF.EXTERNAL_SYS_NUM AS "FED_NM",  UTV80.TINWSF.ST_ASGN_IDENT_CD AS "ST_ID",  UTV80.TINWSF.TYPE_CODE,  UTV80.TINWSYS.NAME AS "SYS_NM", UTV80.TINWSYS.D_PRIN_CNTY_SVD_NM AS "COUNTY", UTV80.TINWSF.NAME AS "FAC_NM",  UTV80.TINWSF.TINWSYS_IS_NUMBER AS "SY_NBR", UTV80.TINLOC.LATITUDE_MEASURE AS "Y",  UTV80.TINLOC.LONGITUDE_MEASURE AS "X",  UTV80.TINLOC.VERTICAL_MEASURE AS "Z", UTV80.TSASAMPL.COLLLECTION_END_DT AS "DTE", UTV80.TSAANLYT.NAME AS "ANLY_NM",  UTV80.TSASAR.CONCENTRATION_MSR AS "CONC_MSR", UTV80.TSASAR.TSASAR_IS_NUMBER AS "ID_NUM",  UTV80.TSASAR.UOM_CODE,   UTV80.TSASAR.DETECTN_LIMIT_NUM AS "DET_LIM", UTV80.TSASAR.DETECTN_LIM_UOM_CD AS "DET_UOM"  FROM  UTV80.TINWSF  INNER JOIN  UTV80.TINWSYS ON UTV80.TINWSF.TINWSYS_IS_NUMBER = UTV80.TINWSYS.TINWSYS_IS_NUMBER  INNER JOIN  UTV80.TINLOC ON UTV80.TINWSF.TINWSF_IS_NUMBER = UTV80.TINLOC.TINWSF_IS_NUMBER  INNER JOIN  UTV80.TSASMPPT ON UTV80.TINWSF.TINWSF_IS_NUMBER = UTV80.TSASMPPT.TINWSF0IS_NUMBER  INNER JOIN  UTV80.TSASAMPL ON UTV80.TSASMPPT.TSASMPPT_IS_NUMBER = UTV80.TSASAMPL.TSASMPPT_IS_NUMBER  INNER JOIN  UTV80.TSASAR ON UTV80.TSASAMPL.TSASAMPL_IS_NUMBER = UTV80.TSASAR.TSASAMPL_IS_NUMBER  INNER JOIN  UTV80.TSAANLYT ON UTV80.TSASAR.TSAANLYT_IS_NUMBER = UTV80.TSAANLYT.TSAANLYT_IS_NUMBER WHERE  UTV80.TINWSYS.D_PRIN_CNTY_SVD_NM LIKE '%CACHE COUNTY%' AND (UTV80.TINWSYS.NAME IN('%PROVID%','%MILL%','%LOG%','%NIB%',   ORDER BY UTV80.TINWSF.ST_ASGN_IDENT_CD
Read in the queried SDWIS data and make a StationId and StationName field.  Make field names consistent with those applied to WQP data above so that compilation is easier later.
End of explanation
"""


print sorted(list(SDWIS.Param.unique()))

def parnormSDWIS(x):
    p = str(x[0]).rstrip().lstrip().lower()
    u = str(x[2]).rstrip().lstrip().lower()
    if p == 'nitrate':
        return 'Nitrate', x[1]*4.427, 'mg/l'
    elif p == 'nitrite':
        return 'Nitrite', x[1]*3.285, 'mg/l'
    elif p == 'nitrogen-ammonia as (n)':
        return 'Ammonium', x[1]*1.288, 'mg/l'
    elif u == 'ug/l':
        return x[0], x[1]/1000, 'mg/l'
    else:
        return x[0], x[1], str(x[2]).rstrip()
    
SDWIS['Param'], SDWIS['ResultValue'], SDWIS['Unit'] = zip(*SDWIS[['Param','ResultValue','Unit']].apply(lambda x: parnormSDWIS(x),1))

"""
Explanation: Normalize units and nutrient data so that they are consistent with the WQP data. This includes standardizing ug/l to mg/l
End of explanation
"""


SDWIS.drop(["FED_NM", "DET_UOM", "UOM_CODE","ANLY_NM", "FAC_NM", "ST_ID", 
            "SYS_NM", "COUNTY", "SY_NBR", "Z", "ID_NUM"],inplace=True, axis=1)

"""
Explanation: Drop unneeded SDWIS fields to save memory and reduce confusion.
End of explanation
"""


SDWISPmatch = {"Ph":"pH","Tds":"TDS","Nitrogen-Ammonia As (N)":"Nitrogen-Ammonia as (N)",
              "Hydroxide As Calcium Carbonate":"Hydroxide as Calcium Carbonate",
               "Bicarbonate As Hco3":"Bicarbonate as HCO3"}
SDWIS["Param"] = SDWIS["Param"].apply(lambda x: SDWISPmatch.get(x,x))
SDWIS["StationName"] = SDWIS["StationName"].apply(lambda x: x.replace("Wtp","WTP"))

SDWIS["ResultValue"] = SDWIS[["ResultValue","Unit"]].apply(lambda x: x[0]/1000 if x[1]=="ug/L" else x[0], 1)

"""
Explanation: Rename chemical parameters in the SDWIS Param field to match those of the WQP data.
End of explanation
"""


SDWISType = {"SP":"Spring","WL":"Well","TP":"Facility Other","IN":"Stream","CC":"Connection","WH":"Well"}
SDWIS.StationType = SDWIS.StationType.apply(lambda x: SDWISType.get(x,x),1)

"""
Explanation: Make station types consistent with the WQP data.
End of explanation
"""


SDWISSta = SDWIS.drop([u'SampleDate', u'ResultValue', u'MDL', u'SampleId', u'Unit', u'MDLUnit', u'Param'], axis=1)
SDWISSta.drop_duplicates(inplace=True)

"""
Explanation: SDWIS facility type code (FacTypeCode): CC, Consecutive_connection;
CH, Common_headers; CS, Cistern; CW, Clear_well;
DS, Distribution_system_zone; IG, Infiltration_gallery; IN,
Intake; NP, Non-piped; OT, Other; PC, Pressure_control;
PF, Pump_facility; RC, Roof_catchment; RS, Reservoir; SI,
Surface_impoundment; SP, Spring; SS, Sampling_station; ST,
Storage; TM, Transmission_main; TP, Treatment_plant; WH,
Well_head; WL, Well.
Create a SDWIS stations file from the SDWIS data. Drop unneeded fields from the station file. Remove duplication stations.
End of explanation
"""


SDWISRes = SDWIS.drop([u'StationType', u'Lat_Y', u'Lon_X', u'StationName', u'Elev'], axis=1)

"""
Explanation: Create SDWIS results file from the SDWIS data. Drop unneeded fields from the results file.  These are fields that are in the station field and apply to stations.
End of explanation
"""


SDWISRes["SampMedia"] = "Groundwater"

"""
Explanation: Create a sample media field and populate it with the value Groundwater.
End of explanation
"""


AGStat = pd.read_csv(AGStationsFile)
AGStat["StationType"] = "Well"
AGStatFields = {"SITEID":"StationId","FINISHEDDE":"Depth","POINT_Y":"Lat_Y",
                "POINT_X":"Lon_X","ELEV_FT":"Elev","ACCURACY":"HorAcc"}
AGStat.rename(columns=AGStatFields,inplace=True)

"""
Explanation: UDAF
UDAF Stations
Import Utah Department of Food and Agriculture data from the data file. These data were compiled from <a href=http://ag.utah.gov/conservation-environmental/ground-water.html>reports available of the UDAF website</a>. Once the data are imported, rename the fields to match the above SDWIS and WQP data.
End of explanation
"""


AGStat.drop(["OBJECTID_1", "OBJECTID", "PUB_YR", "SAMPLENO", "WLDATE", "WLDEPTH"], inplace=True, axis=1)

"""
Explanation: Drop unneeded fields to save memory.
End of explanation
"""


AGStat.StationId = AGStat.StationId.apply(lambda x: "UDAF-"+str(int(x)).zfill(5),1)

"""
Explanation: Add UDAF prefix to the station identification field (StationId) to make station ids unique.
End of explanation
"""


names = ["SampleId","ResultValue", "ParAbb", "Unit", "Param", "MDL","BelowLim","TestNo",
         "StationId","SampleDate","SampYear"]
AGRes = pd.read_csv(AGResultsFile, names=names, index_col=10)

"""
Explanation: UDAF Results
Import data Utah Department of Food and Agriculture data from the data file. These data were compiled from reports available of the UDAF website. Once the data are imported, rename the fields to match the above SDWIS and WQP data.
End of explanation
"""


AGRes["DetectCond"] = AGRes["BelowLim"].apply(lambda x: 'Not Detected' if x=='Y' else np.nan,1)

"""
Explanation: Create a detection condition field and populate it based on values in the imported data.
End of explanation
"""


AGRes.ResultValue = AGRes[["BelowLim","ResultValue"]].apply(lambda x: np.nan if x[0]=="Y" or x[1] == 0.0 else x[1], 1)

"""
Explanation: Fill null result values with zeros when data are reported as below detection limit.
End of explanation
"""


def parnormAG(x):
    p = str(x[0]).rstrip().lstrip().lower()
    u = str(x[2]).rstrip().lstrip().lower()
    if p == 'nitrate-n':
        return 'Nitrate', x[1]*4.427, 'mg/l'
    elif u == 'ug/l':
        return x[0], x[1]/1000, 'mg/l'
    else:
        return x[0], x[1], str(x[2]).rstrip()
    
AGRes['Param'], AGRes['ResultValue'], AGRes['Unit'] = zip(*AGRes[['Param','ResultValue','Unit']].apply(lambda x: parnormAG(x),1))

AGRes.Unit.unique()

AGRes.dropna(subset=["StationId","ResultValue"], how="any", inplace=True)

AGRes.StationId = AGRes.StationId.apply(lambda x: "UDAF-"+str(int(x)).zfill(5),1)

AGStAv = list(AGStat.StationId.values)
AGRes = AGRes[AGRes.StationId.isin(AGStAv)]
AGRes["SampMedia"] = "Groundwater"

AGStat['OrgId']='UDAF'

"""
Explanation: Make data consistent by cleaning up parameter descriptions.
End of explanation
"""


STORLegSta = pd.read_table(STORLegStatFile, skiprows=[1])
STORLegRes = pd.read_table(STORLegResFile, skiprows=[1])
STORParam = pd.read_table(STORParamFile)

"""
Explanation: STORET Legacy
Legacy EPA data are kept in the <a href=ftp://ftp.epa.gov/storet/exports/>STORET Legacy Database</a>.
End of explanation
"""


rescol = list(STORLegRes.columns)
j = []
for i in rescol:
    j.append(i.rstrip("\t").rstrip().lstrip().replace(" ",""))
resdict = dict(zip(rescol,j))
STORLegRes.rename(columns=resdict,inplace=True)

statcol = list(STORLegSta.columns)
k = []
for i in statcol:
    k.append(i.rstrip("\t").rstrip().lstrip().replace(" ",""))
statdict = dict(zip(statcol,k))
STORLegSta.rename(columns=statdict,inplace=True)

STORLegRes["SampleDate"] = STORLegRes[["StartDate","StartTime"]].apply(lambda x: datetimefix(x,"%Y-%m-%d %H:%M"),1)
STORLegRes = STORLegRes[STORLegRes.SecondaryActivityCategory.isin(['Water',np.nan])]
STORParamDict = dict(zip(STORParam['Parameter No.'].values, STORParam['Full Name'].values))
STORLegRes.Param = STORLegRes.Param.apply(lambda x: STORParamDict.get(x),1)

STORResField = {"Agency":"OrgId","AgencyName":"OrgName","Station":"StationId","SampleDepth":"SampDepth"}
STORLegRes.rename(columns=STORResField,inplace=True)

STORLegRes.drop(["StateName", "CountyName", "HUC", "EndDate", "UMK", "CS", "ReplicateNumber",
                 "COMPOSITE_GRAB_NUMBER","CM","PrimaryActivityCategory","PrimaryActivityCategory",
                 "SecondaryActivityCategory",
                 "EndTime", "StartDate", "StartTime", "Latitude", "Longitude"],inplace=True,axis=1)

STORLegRes["SampleId"] = STORLegRes[["StationId","SampleDate"]].apply(lambda x: str(x[0]) + "-" + str(x[1]),1 )
STORLegRes["StationId"] = STORLegRes["StationId"].apply(lambda x: "EPALeg-" + x, 1)

STORLegRes.Param = STORLegRes.Param.apply(lambda x: str(x).title(),1)

STORLegRes.columns

def parnormSTOR(x):
    p = str(x[0]).rstrip().lstrip().lower()
    if p ==  'nitrate nitrogen, total (mg/L as n)' or p== 'nitrate nitrogen, total':
        return 'Nitrate', x[1]*4.427, 'mg/l'
    elif p == 'nitrite nitrogen, total (mg/l as n)':
        return 'Nitrite', x[1]*3.285, 'mg/l'
    elif p == 'nitrogen, ammonia, total (mg/l as n)':
        return 'Ammonium', x[1]*1.288, 'mg/l'
    elif p == 'sulfate (as s) whole     water, mg/L':
        return 'Sulfate', x[1]*2.996, 'mg/l'
    elif p in ('phosphorus, dissolved orthophosphate (mg/l as p)'):
        return 'Phosphate', x[1]*3.066, 'mg/l'
    else:
        return x[0], x[1], np.nan

STORLegRes['Param'], STORLegRes['ResultValue'], STORLegRes['Unit'] = zip(*STORLegRes[['Param','ResultValue']].apply(lambda x: parnormSTOR(x),1))

STORKeepers = ['Temperature, Water (Degrees Centigrade)',
               'Temperature, Water (Degrees Fahrenheit)', 
               'Specific Conductance,Field (Umhos/Cm @ 25C)',
               'Specific Conductance (Umhos/Cm @ 25C)', 
               'Sulfate (As S) Whole     Water, Mg/L',
               'Oxygen, Dissolved                             Mg/L',
               'Oxygen, Dissolved, Percent Of Saturation         %',
               'Bod, 5 Day, 20 Deg C                          Mg/L',
               'Ph (Standard Units)', 'Ph, Lab, Standard Units                         Su',
               'Carbon Dioxide (Mg/L As Co2)', 'Alkalinity,Total,Low Level Gran Analysis     Ueq/L',
               'Alkalinity, Total (Mg/L As Caco3)', 'Bicarbonate Ion (Mg/L As Hco3)', 'Carbonate Ion (Mg/L As Co3)',
               'Nitrogen, Ammonia, Total (Mg/L As N)', 'Ammonia, Unionzed                      (Mg/L As N)',
               'Nitrite Nitrogen, Total (Mg/L As N)', 'Ammonia, Unionized (Calc Fr Temp-Ph-Nh4)  (Mg/L)',
               'Nitrate Nitrogen, Total (Mg/L As N)', 'Nitrogen, Kjeldahl, Total, (Mg/L As N)',
               'Nitrite Plus Nitrate, Total 1 Det. (Mg/L As N)', 'Phosphorus (P), Water, Total Recoverable      Ug/L',
               'Phosphorus, Total (Mg/L As P)', 'Phosphorus, Dissolved Orthophosphate (Mg/L As P)',
               'Carbon, Dissolved Organic (Mg/L As C)',
               'Carbon, Dissolved Inorganic (Mg/L As C)',
               'Hardness, Total (Mg/L As Caco3)', 'Calcium (Mg/L As Caco3)',
               'Calcium, Dissolved (Mg/L As Ca)',
               'Magnesium, Dissolved (Mg/L As Mg)',
               'Sodium, Dissolved (Mg/L As Na)',
               'Potassium, Dissolved (Mg/L As K)',
               'Chloride, Dissolved In Water           Mg/L',
               'Sulfate, Dissolved (Mg/L As So4)',
               'Fluoride, Dissolved (Mg/L As F)',
               'Silica, Dissolved (Mg/L As Si02)',
               'Arsenic, Dissolved  (Ug/L As As)', 'Arsenic, Total (Ug/L As As)',
               'Barium, Dissolved (Ug/L As Ba)', 'Barium, Total (Ug/L As Ba)',
               'Beryllium, Total (Ug/L As Be)', 'Boron, Dissolved (Ug/L As B)',
               'Boron, Total (Ug/L As B)', 'Cadmium, Dissolved (Ug/L As Cd)',
               'Cadmium, Total (Ug/L As Cd)', 'Chromium, Dissolved (Ug/L As Cr)',
               'Chromium, Hexavalent (Ug/L As Cr)', 'Chromium, Total (Ug/L As Cr)',
               'Copper, Dissolved (Ug/L As Cu)', 'Copper, Total (Ug/L As Cu)',
               'Iron, Dissolved (Ug/L As Fe)', 'Lead, Dissolved (Ug/L As Pb)',
               'Lead, Total (Ug/L As Pb)', 'Manganese, Total (Ug/L As Mn)',
               'Manganese, Dissolved (Ug/L As Mn)', 'Thallium, Total (Ug/L As Tl)',
               'Nickel, Dissolved (Ug/L As Ni)', 'Nickel, Total (Ug/L As Ni)',
               'Silver, Dissolved (Ug/L As Ag)', 'Silver, Total (Ug/L As Ag)',
               'Zinc, Dissolved (Ug/L As Zn)', 'Zinc, Total (Ug/L As Zn)',
               'Antimony, Total (Ug/L As Sb)', 'Aluminum, Total (Ug/L As Al)',
               'Selenium, Dissolved (Ug/L As Se)', 'Selenium, Total (Ug/L As Se)',
               'Tritium (1H3),Total (Picocuries/Liter)',
               'Hardness, Ca Mg Calculated (Mg/L As Caco3)',
               'Chlorine, Total Residual (Mg/L)',
               'Residue,Total Filtrable (Dried At 180C),Mg/L',
               'Nitrate Nitrogen, Dissolved (Mg/L As No3)', 'Iron (Ug/L As Fe)',
               'Phosphorus, Total, As Po4 - Mg/L', 'Mercury, Total  (Ug/L As Hg)']
STORLegRes = STORLegRes[STORLegRes.Param.isin(STORKeepers)]

def parsplit(x,p):
    x = str(x).rstrip().lstrip()
    if p == "Un":
        z = -1
        x = str(x).replace("Mg/L", "mg/l")
        x = str(x).replace("Ug/L", "ug/l")
        x = str(x).replace("o", "O")
        x = str(x).replace("c", "C")
        x = str(x).replace("TOtal ReCOverable","Total Recoverable")        
        x = str(x).replace("UmhOs", "umhos")
        x = str(x).replace("TOtal","Total")
    elif p== "Par":
        z = 0
        x = str(x).replace(", Standard Units","")
        x = str(x).replace(", Unionized","")
        x = str(x).replace(", Unionzed","")
        x = str(x).replace(",Low Level Gran Analysis","")
        x = str(x).replace(" Ion","")
        x = str(x).replace(",Total",", Total")
        if x == "Ph" or x == "Ph, Lab":
            x = str(x).replace("Ph","pH")
    if "(" in x:
        x = str(x).replace(" As ", " as ")
        return str(x).split(" (")[z].rstrip(")").rstrip().lstrip()
    else:
        return str(x).split("  ")[z].rstrip().lstrip()

def splitmore(x):
    if "NO3" in x:
        return x
    elif " as " in x:
        return x.split(" as ")[0]
    elif x == "As S) WhOle     Water, mg/l" or x == "Dried At 180C),mg/l" or x=="PhOsphOrus, Total, As PO4 - mg/l":
        return "mg/l"
    elif x == "P), Water, Total Recoverable      ug/l":
        return "ug/l"
    else:
        return x

def unitconv(x):
    if x[1]=="ug/l":
        return x[0]/1000
    elif x[1]=="Degrees Fahrenheit":
        return (float(x[0])-32.0)*(5.0/9.0)
    else:
        return x[0]

STORLegRes["Unit"] = STORLegRes["Param"].apply(lambda x: parsplit(x,"Un"), 1)
STORLegRes["Param"] = STORLegRes["Param"].apply(lambda x: parsplit(x,"Par"), 1)
STORLegRes["Unit"] = STORLegRes["Unit"].apply(lambda x: splitmore(x), 1)
STORLegRes["ResultValue"] = STORLegRes[["ResultValue","Unit"]].apply(lambda x: unitconv(x), 1)
STORLegRes["Unit"] = STORLegRes["Unit"].apply(lambda x: "mg/l" if x=="ug/l" else x, 1)
STORLegRes["Unit"] = STORLegRes["Unit"].apply(lambda x: "Degrees Centigrade" if x=="Degrees Fahrenheit" else x, 1)

STORStaField = {"Agency":"OrgId","AgencyName":"OrgName","Station":"StationId", "DepthUnits":"DepthUnit",
                "Latitude":"Lat_Y", "Longitude":"Lon_X", "HUC":"HUC8", "StationDepth":"Depth"}
STORLegSta.rename(columns=STORStaField,inplace=True)

STORLegSta.columns

STORLegSta.drop(["RchmileSegment", "MilesUpReach", "Rchonoff", "Description", "G", "S", "StationAlias",
                 "Rchname", "StateName", "CountyName"], inplace=True, axis=1)

STORLegSta.StationType = STORLegSta.StationType.apply(lambda x: str(x).rstrip(" ").strip("/SUPPLY").split("/")[-1].title(),1)

LegTypeDict = {"We":"Well"}
STORLegSta.StationType = STORLegSta.StationType.apply(lambda x: LegTypeDict.get(x,x),1)

STORLegSta.StationId = STORLegSta["StationId"].apply(lambda x: "EPALeg-" + x, 1)

"""
Explanation: Parse choppy text data from the STORET Legacy database.
End of explanation
"""


UGSfield = pd.read_excel(fielddata,"FieldChem") #Field data
UGSNO3 = pd.read_excel(fielddata,"Nitrate") #Nitrate data provided by Millville City

UGS = pd.read_csv(UGSFile, engine="python")
UGS["StationId"] = UGS["SITE"].apply(lambda x:"UGS-"+str(x).zfill(4),1)

UGSSta = UGS.drop([u'OBJECTID_1',u'SITE', u'TDS', u'Temp', u'Cond', u'CO2', u'HCO3', 
                  u'CO3',u'Na', u'pH', u'Ca', u'SO4', u'NO3', u'As_', u'Cl', u'K',
                  u'Mg', u'Hard', u'NH4'], axis=1)

UGSRe = UGS.drop([u'OBJECTID_1',u'SITE',u'StationType', u'Geology', u'Elev', u'Lat_Y', u'Lon_X', u'StationName', 
                  u'OrgId', u'WRNUM', u'SITE', u'UTM_X', u'UTM_Y', u'Depth_ft'], axis=1)

UGSRe["SampleId"] = UGSRe.index

UGSRe.reset_index(inplace=True)
UGSRe.set_index(["StationId","SampleId"], inplace=True)

UGSRe.drop(UGSRe.columns[0],inplace=True,axis=1)

UGSStack = UGSRe.stack().to_frame()
UGSStack.columns = ["ResultValue"]

UGSStack.reset_index(inplace=True)

UGSStack.columns=["StationId","SampleId","Param","ResultValue"]

def unitcon(x):
    if x=="pH":
        return ""
    elif x=="Temp":
        return "C"
    elif x=="Cond":
        return "uS/cm"
    else:
        return "mg/l"
    
UGSStack["Unit"] = UGSStack["Param"].apply(lambda x: unitcon(x),1)
UGSStack["ParAbb"] = UGSStack["Param"]
UGSStack["OrgId"] = "UGS"
UGSStack["OrgName"] = "Utah Geological Survey"
UGSStack["ResultValue"] = UGSStack[['Param','ResultValue']].apply(lambda x: x[1]*1.288 if x[0]=='Ammonia as N' else x[1],1)
UGSStack["Param"] = UGSStack['Param'].apply(lambda x: 'Ammonia' if x=='Ammonia as N' else x, 1)
UGSStack["ResultValue"] = UGSStack[['Param','ResultValue']].apply(lambda x: x[1]*3.066 if x[0]=='Phosphate, Tot. Dig. (as P)' else x[1],1)
UGSStack["Param"] = UGSStack['Param'].apply(lambda x: 'Phosphate' if x=='Phosphate, Tot. Dig. (as P)' else x, 1)

"""
Explanation: UGS Data
End of explanation
"""


SLSampMatch = pd.read_excel(fielddata,"StateLabMatch")
SLStat = pd.read_excel(fielddata,"Stations")
#SLStat = pd.merge(SLSampMatch, SLStations, on='StationId', how='outer')
#SLStat.reset_index(inplace=True)

SLStat

SL0 = pd.read_table(statelabresults0, sep="\t", lineterminator="\n", error_bad_lines=False)
SL0 = SL0[SL0['Collector']=='PI']
SL1 = pd.read_table(statelabresults1, sep="\t", lineterminator="\n", error_bad_lines=False)
SL1 = SL1[SL1['Collector']=='PI']
SL2 = pd.read_table(statelabresults2, sep="\t", lineterminator="\n", error_bad_lines=False)
SL2 = SL2[SL2['Collector']=='PI']

SL = pd.concat([SL0,SL1,SL2])

SL["OrgId"] = "UGS"
SL["OrgName"] = "Utah Geological Survey"
SL['DetectCond'] = SL['Problem#Identifier'].apply(lambda x: 'Not Detected' if str(x).rstrip()=='<' else np.nan,1)
SL['SampleDate'] = SL[['Sample#Date','Sample#Time']].apply(lambda x: datetimefix(x,"%m/%d/%y %H:%M"),1)

SLHead = {'Sample#Number':'SampleId', 'Param#Description':'Param', 'Result#Value':'ResultValue','Units':'Unit', 
          'Lower#Report#Limit':'MDL','Method#ID':'SampMeth','Analysis#Date':'AnalysisDate'}
SL.rename(columns=SLHead,inplace=True)

SL['Sample#Description'].unique()

SL.drop([u'Lab#Code', u'Station#ID', u'Source#Code', u'Sample#Date',
       u'Sample#Time', u'Sample#Type', u'Cost#Code', u'Billing#Code',
       u'Agency#Bill#Code', u'Trip#ID', u'Sample#Description', u'Collector',
       u'Sample#Recieved#Date', u'Chain#of#Custody#Ind.', u'Replicate#Number',
       u'Sample#Comment', u'Method#Number', u'Method#Agency',
       u'Method#Description', u'Param#Number', u'CAS#Number',
       u'Matrix#Number', u'Matrix#Description', u'Preparation#Date',
       u'Problem#Identifier', u'Result#Code', 
       u'Upper#Quant#Limit', u'Method#Detect#Limit',
       u'Confidence#Limit', u'%#Confidence#Limit',u'Dilution#Factor',
       u'Batch#Number',u'Comment#Number', u'Comment#Text'], inplace=True, axis=1)

SL.columns

SLRes = pd.merge(SL, SLSampMatch, on='SampleId', how='left')

SLStat.drop_duplicates(inplace=True)

def SLparnorm(x):
    p = str(x[0]).rstrip().lstrip().lower()
    u = str(x[2]).rstrip().lstrip().lower()
    if p == 'nitrate nitrogen, total (mg/l as n)':
        return 'Nitrate', x[1]*4.427, 'mg/l'
    elif p == 'nitrite nitrogen, total (mg/l as n)':
        return 'Nitrite', x[1]*3.285, 'mg/l'
    elif p == 'ammonia as n':
        return 'Ammonium', x[1]*1.288, 'mg/l'
    elif p == 'sulfate (as s) whole     water, mg/L':
        return 'Sulfate', x[1]*2.996, 'mg/l'
    elif p in ('phosphate, tot. dig. (as p)', 'phosphate-phosphorus as p','orthophosphate as p'):
        return 'Phosphate', x[1]*3.066, 'mg/l'
    elif u == 'ug/l':
        return x[0], x[1]/1000, 'mg/l'
    else:
        return x[0], x[1], str(x[2]).rstrip()

def MDLfix(x):
    u = str(x[1]).rstrip().lstrip().lower()
    if np.isfinite(x[2]):
        return x[0]
    elif u=='ug/l':
        return x[0]/1000
    else:
        return x[0]

SLRes['MDL'] = SLRes[['MDL','Unit','ResultValue']].apply(lambda x: MDLfix(x),1)
SLRes['Param'], SLRes['ResultValue'], SLRes['Unit'] = zip(*SLRes[['Param','ResultValue','Unit']].apply(lambda x: parnorm(x),1))

SLRes.StationId.unique()

"""
Explanation: State Lab
These are raw data results sent to the UGS via Tab-delimited tables from the Utah State Health Laboratory.  They make up the bulk of results of data collected for this study. They are supplimented with field data translated to spreadsheets.
End of explanation
"""


Res = pd.concat([STORLegRes,AGRes,SDWISRes,WQP,UGSStack,SLRes,UGSfield,UGSNO3])

Res = Res[~Res["Unit"].isin(['ueq/L','Ueq/L','ueq/l','tons/ac ft','tons/day','meq/L'])]
Res = Res[~Res["Param"].isin(["Heptachlorobiphenyl", "Hydrocarbons", "Hydroxide", "Ionic strength",
                              "Floating debris, severity", "Carbon Tetrachloride", "Trichlorobiphenyl",
                              "Vinyl Chloride", "True color", "Color", "Trash, Debris, Floatables",
                              "Total volatile solids", "Temperature, air", "Residue, Total Filtrable",
                              "Pentachlorobiphenyl", "Odor threshold number", "Odor, atmospheric",
                              "Instream features, est. stream width", "Hydroxide", 
                              "Light, transmissivity","Algae, floating mats (severity)"])]
len(Res)

Res[["Param","Unit","USGSPCode"]].drop_duplicates(subset=["Param","Unit"]).sort_values(by=["Param"]).to_clipboard()

Stat = pd.concat([STORLegSta, AGStat, SDWISSta, WQPStat, SLStat, UGSSta])

parmatch = pd.read_excel(rootname + "Aquachem.xlsx")

parmatchdict = dict(zip(parmatch.Param.values, parmatch.ParrAbb.values))
Res["ParAbb"] = Res[["ParAbb","Param"]].apply(lambda x: parmatchdict.get(x[1],x[0]),1)

results = Res.dropna(subset=["StationId","Param","SampleId"], how="any")

Stat.loc[:,"StationName"] = Stat["StationName"].apply(lambda x: str(x).strip().lstrip().rstrip(),1)
Stat.loc[:,"StationId"] = Stat["StationId"].apply(lambda x: str(x).strip().lstrip().rstrip(),1)
Res.loc[:,"StationId"] = Res["StationId"].apply(lambda x: str(x).strip().lstrip().rstrip(),1)

results.loc[:,"Unit"] = results[["ParAbb","Unit"]].apply(lambda x: "C" if x[0]=="Temp" else x[1],1)
results.loc[:,"Unit"] = results[["ParAbb","Unit"]].apply(lambda x: "umhos/cm" if x[0]=="Cond" else x[1],1)
results.loc[:,"Unit"] = results[["ParAbb","Unit"]].apply(lambda x: "" if x[0]=="pH" else x[1],1)

results.drop(["AnalysisDate","AnalytMeth","SampType","AnalytMethId", "BelowLim", "StationName",
              "MethodDescript", "LabComments", "LabName", "LimitType", "ProjectId", "QualCode",
              "OrgName","R", "ResultComment","ResultStatus","SampComment", "SampEquip", 
              "SampDepthRef", "SampDepthU","SampDepth","SampType", "USGSPCode",
              "SampMeth", "SampMethName","SampYear","TestNo"],inplace=True,axis=1)

"""
Explanation: Combine Data
End of explanation
"""


NDs = {'Not Detected':'<', 'Present Above Quantification Limit':'<', 'ND             ':'<', '*Present >QL   ':'>',
       'Present Below Quantification Limit':'<', '*Non-detect    ':'<', 'Detected Not Quantified':'<', 
       'Systematic Contamination':'<'}
results.DetectCond = results.DetectCond.apply(lambda x: NDs.get(x,np.nan),1)

def is_nan(x): 
    '''
    this function identifies nan values
    Source: http://stackoverflow.com/questions/944700/how-to-check-for-nan-in-python
    '''
    try: 
        return math.isnan(x) 
    except: 
        return False 

def detected(x):
    '''
    Finds nondetects and fixes units and values
    '''
    if x[1]=='<' and np.isfinite(x[0]):
        return x[1]+str(x[0])       
    elif x[1]=='<' and np.isfinite(x[2]):
        if str(x[3]).rstrip().lower() == 'ug/l':
            return x[1]+str(x[2]/1000)
        else:
            return x[1]+str(x[2])
    else:
        return x[0]

results.ResultValue = results[['ResultValue','DetectCond','MDL','MDLUnit']].apply(lambda x: detected(x),1)

def MDLfill(x):
    if x[0] <= 0 and x[1]>0:
        return 0
    elif x[2] == '<':
        return 0
    elif x[0] < x[1]:
        return 0
    else:
        return 1
results.loc[:,'ResValue'] = pd.to_numeric(results['ResultValue'], errors='coerce')    
results.loc[:,'Censored'] = results[['ResValue','MDL','DetectCond']].apply(lambda x: MDLfill(x),1)

matchDict = {'414143111495501':'USGS-414143111495501','414115111490301':'USGS-414115111490301',
             'SDWIS3117.0WS004':'USGS-414115111490301', 
             'EPALeg-0301203':'414029111483501','SDWIS3116.0WS003':'414029111483501',
             'EPALeg-0301201':'USGS-414024111481101','SDWIS5435.0WS001':'USGS-414024111481101',
             '414024111481101':'USGS-414024111481101','EPALeg-0300101':'SDWIS5411.0WS001',
             'EPALeg-0300102':'SDWIS5412.0WS002', 'EPALeg-0300103':'SDWIS5413.0WS003', 
             'UGS-107.5':'SDWIS3143.0WS001','UDAF-01492':'UGS-0412', 'UDAF-03165':'UGS-106.5',
             'SDWIS3126.0WS002':'USGS-414216111485201', 'EPALeg-0301702':'USGS-414216111485201',
             'EPALeg-0301901':'USGS-414328111493001', 'SDWIS3131.0WS001':'USGS-414328111493001', 
             'EPALeg-0301005':'SDWIS3112.0WS005', 'EPALeg-0301002':'USGS-414417111484301', 
             'SDWIS3109.0WS002':'USGS-414417111484301', 'SDWIS3113.0WS006':'USGS-414459111493601', 
             'SDWIS3127.0WS003':'414213111493101', 'SDWIS3159.0WS003':'SDWIS3157.0WS001','UDAF-01500':'UGS-63.5',
             'SDWIS3111.0WS004':'USGS-414441111490701', 'EPALeg-0301904':'SDWIS3133.0WS004',
            'EPALeg-0301004':'USGS-414441111490701','EPALeg-0301502':'SDWIS3118.0WS002',
             'UDAF-01589':'UDAF-01568','UDAF-01586':'UDAF-01566','UGS-0050':'UDAF-01566',
            'EPALeg-0300104':'SDWIS3088.0WS004', 'UDAF-01585':'UGS-0032', 'UDAF-01565':'UGS-0032',
            'EPALeg-0300201':'SDWIS3091.0WS001', 'EPALeg-0300204':'SDWIS3094.0WS004',
            'EPALeg-0301803':'SDWIS3129.0WS003','EPALeg-0300405':'SDWIS5418.0WS005',
            'EPALeg-0300404':'SDWIS5417.0WS004', 'EPALeg-0300403':'SDWIS5416.0WS003',
            'SDWIS5439.0WS003':'SDWIS5416.0WS003', 'SDWIS5460.0WS003':'SDWIS5416.0WS003',
            'SDWIS5414.0WS001':'SDWIS5458.0WS001', 'SDWIS5437.0WS001':'SDWIS5458.0WS001',
            'EPALeg-0308601':'USGS-415828111460001', 'SDWIS5487.0WS001':'USGS-415828111460001',
            'SDWIS5430.0WS002':'USGS-415828111460001', 'SDWIS5423.0WS003':'USGS-415828111460001',
            'SDWIS5421.0WS001':'USGS-415836111464701', 'EPALeg-0304901':'SDWIS5479.0WS001', 
            'EPALeg-0303201':'SDWIS5470.0WS001', 'SDWIS5432.0WS001':'USGS-414535111423001',
            'EPALeg-0303001':'SDWIS5469.0WS001','EPALeg-0307701':'SDWIS5485.0WS001',
            'EPALeg-0308301':'SDWIS5486.0WS001', 'EPALeg-0301501':'SDWIS5445.0WS001',
            'EPALeg-0300701':'USGS-415120111440001', 'SDWIS5424.0WS001':'USGS-415120111440001',
            'EPALeg-0302001':'SDWIS5455.0WS001', 'EPALeg-0301101':'SDWIS5433.0WS001', 
            'EPALeg-0301102':'SDWIS5434.0WS002'}


Stat.loc[:,'StationId'] = Stat['StationId'].apply(lambda x: matchDict.get(x,x),1)
results.loc[:,'StationId'] = results['StationId'].apply(lambda x: matchDict.get(x,x),1)             
results.loc[:,'SampleDate'] = pd.to_datetime(results.SampleDate)

def depthFill(x):
    if x > 0:
        return x
    
def depthUnitFill(x):
    if x > 0:
        return 'ft'

Stat.Depth = Stat['Depth_ft'].apply(lambda x: depthFill(x),1)
Stat.DepthUnit = Stat['Depth_ft'].apply(lambda x: depthUnitFill(x),1)

WINdict = {'SDWIS3180.0WS001':435116, 'UGS-47.5':32700, 'UDAF-01566':30211, 'UGS-46.5':12420, 
           'USGS-414525111503705':427268, 'UDAF-01569':28327, 'SDWIS3112.0WS005':2694,
           'UDAF-03162':434818, 'USGS-414328111493001':2823, 'USGS-414332111491001':2836,
           'SDWIS3133.0WS004':2848, 'UGS-91.5':28647, 'UGS-95.5':35814, 'SDWIS3128.0WS004':18590,
           'UGS-0102':426853, 'USGS-414115111490301':2722, 'UT4140521114843201':32975, '414029111483501':2721,
           'SDWIS3088.0WS004':2741,'UGS-63.5':7126, 'UGS-0084':9639, 'USGS-414134111544701':434098, 'UGS-0070':35061,
           'UGS-0029':32851,'UGS-0030':26663,'UGS-0034':29110, 'UGS-0055':3728, 'UGS-61.5':9280, 
           'SDWIS3129.0WS003':2816, 'UGS-0043':29329, 'UGS-0889':24493, 'UGS-44.5':28333}
 
Stat.WIN = Stat['StationId'].apply(lambda x: WINdict.get(x,x),1)

results.SampleID = results.SampleId.apply(lambda x: str(x).replace(' ',''),1)
results.StationId = results.StationId.apply(lambda x: str(x).replace(' ',''),1)
Stat.StationId = Stat.StationId.apply(lambda x: str(x).replace(' ',''),1)

results.drop_duplicates(subset = ['SampleId','ParAbb'],inplace=True)
Stat.drop_duplicates(subset = ['StationId'],inplace=True)

resultsNoND = results[(~results['DetectCond'].isin(['<','>']))]

"""
Explanation: Clean Up Non Detects
End of explanation
"""


datap = resultsNoND.pivot(index='SampleId', columns='ParAbb', values='ResValue')

datap.dropna(subset=['SO4','Cond','Temp','TDS','pH_field'],how='all',inplace=True)
datap.drop(datap.columns[[0]], axis=1, inplace=True)

results.columns

resdrop = [ 'DetectCond', u'Comment#Number.1',   u'Comment#Text.1', 'ResultValue', 'ResValue',
             'MDL', 'MDLUnit', 'OrgId', 'Param', 'ResultValue', 'SampFrac',
             'SampMedia', 'Unit', 'ParAbb']
resPivot = results.drop(resdrop, axis=1)

datapiv = pd.merge(datap, resPivot, left_index=True, right_on='SampleId',how='left')
datapiv.drop_duplicates(subset=['SampleId'],inplace=True)

"""
Explanation: Pivot Data
End of explanation
"""


def projy(x):
    inProj = Proj(init='epsg:4326') #WGS84
    outProj = Proj(init='epsg:2152') #NAD83(CSRS98) / UTM zone 12N
    x2,y2 = transform(inProj,outProj,x[0],x[1])
    return y2

def projx(x):
    inProj = Proj(init='epsg:4326') #WGS84
    outProj = Proj(init='epsg:2152') #NAD83(CSRS98) / UTM zone 12N
    x2,y2 = transform(inProj,outProj,x[0],x[1])
    return x2

def getelev(x):
    elev = "http://ned.usgs.gov/epqs/pqs.php?x="+str(x[0])+"&y="+str(x[1])+"&units=Meters&output=xml"
    response = urllib2.urlopen(elev)
    html = response.read()
    d = xmltodict.parse(html)
    return float(d['USGS_Elevation_Point_Query_Service']['Elevation_Query']['Elevation'])

Stat.loc[:,'UTM_X'] = Stat[['Lon_X','Lat_Y']].apply(lambda x: projx(x),1)
Stat.loc[:,'UTM_Y'] = Stat[['Lon_X','Lat_Y']].apply(lambda x: projy(x),1)
Stat.loc[:,'Elev'] = Stat[['Lon_X','Lat_Y']].apply(lambda x: getelev(x),1)

pivStats = Stat.drop(['Aquifer', 'ConstDate', 'Depth', 'DepthUnit','AquiferType', 'HorCollMeth', 'Geology',
                         'HoleDUnit', 'HoleDepth', 'HUC8', 'HorAccUnit', 'HoleDUnit', 'SCREENDEPT',
                         'ElevUnit', 'ElevRef', 'ElevAcc', 'ElevMeth','CountyCode', 'ElevAccUnit',
                         'HorAcc', 'StateCode', 'HorRef',
                         'OrgId', 'StationComment'], axis=1)

pivStats.reset_index(inplace=True)
pivStats.set_index("StationId",inplace=True)

pivdata = pd.merge(datapiv, pivStats, left_on="StationId", right_index=True, how='left')

pivdata.drop_duplicates(subset=['SampleId'],inplace=True)

"""
Explanation: Add GIS Information
End of explanation
"""


alkmatch = pivdata[(pivdata['Meas_Alk']>0)&(pivdata['HCO3']>0)]

x = [np.float64(i) for i in alkmatch['Meas_Alk'].values]
y = [np.float64(i) for i in alkmatch['HCO3'].values]

X = sm.add_constant(x)
res = sm.RLM(y,X).fit()
b = res.params[0]
m = res.params[1]
print m
print b

plt.figure()
plt.scatter(x,y)
plt.plot(x, res.fittedvalues, color='red')

def HCO3fix(x):
    if x[0]>0:
        return x[0]
    elif x[1]>0:
        return x[1]*m+b
    else:
        pass

pivdata['HCO3'] = pivdata[['HCO3','Meas_Alk']].apply(lambda x: HCO3fix(x),1)

parlist = ['Ca','Mg','Na','K','Cl','HCO3','CO3','SO4','NO3','NO2','CO2','TDS','Si','Zn_tot','As_tot']

def removeInf(x):
    if x <= 0:
        return np.nan
    else:
        return np.log(x)

for i in parlist:
    if i in pivdata.columns:
        pivdata[i+'Ln'] = pivdata[i].apply(lambda x: removeInf(x),1)

d = {'Ca':0.04990269, 'Mg':0.082287595, 'Na':0.043497608, 'K':0.02557656, 'Cl':0.028206596, 'HCO3':0.016388838, 'CO3':0.033328223, 'SO4':0.020833333, 'NO2':0.021736513, 'NO3':0.016129032}
chemlist = ['Ca','Mg','Na','K','Cl','HCO3','CO3','SO4','NO3','NO2']

for i in chemlist:
    if i in pivdata.columns:
        pivdata[i+'Meq'] = pivdata.loc[:,i] * d[i]

pivdata.drop_duplicates(subset = ['StationId','SampleDate'], inplace=True)

def sumIons(x):
    b = 0
    for i in x:
        if i>0:
            b = i + b
        else:
            b = b
    return b

pivdata['Anions'] = pivdata[['ClMeq','HCO3Meq','SO4Meq','CO3Meq']].apply(lambda x: sumIons(x),1)
pivdata['Cations'] = pivdata[['KMeq','MgMeq','NaMeq','CaMeq']].apply(lambda x: sumIons(x),1)
pivdata['EC'] = pivdata['Anions'] - pivdata['Cations'] 
pivdata['CBE'] = ((pivdata['Cations']-np.abs(pivdata['Anions']))/(pivdata['Cations']+np.abs(pivdata['Anions'])))*100

"""
Explanation: Convert and Balance Samples
End of explanation
"""


#piperdata = pivdata.dropna(subset = ['Ca','Na','Cl','Mg','SO4','HCO3'], how='any')
#piperdata.drop_duplicates(subset=['SampleId'], inplace=True)

print(len(pivdata))

pivgrps = pivdata.groupby(['StationId']).median()
pivGoodData = pivdata[abs(pivdata.CBE)<=5]
pipergrps = pivGoodData.groupby(['StationId']).median()
pipergrps['sampCount'] = pivGoodData.groupby(['StationId'])['CBE'].agg({'cnt':(lambda x: np.count_nonzero(~np.isnan(x)))}).reset_index#squeeze=True

pivgrp = pd.merge(pivgrps, pivStats, left_index=True, right_index=True, how='left')
pipergrp = pd.merge(pipergrps, pivStats, left_index=True, right_index=True, how='left')
pipergrp.drop_duplicates(inplace=True)
pivgrp = pivgrp.reset_index().drop_duplicates(subset=['StationId']).set_index('StationId')

princpiv = pivGoodData[(pivGoodData.SampleDate < datetime.datetime(2014,3,10))&(pivGoodData.UTM_X < 435000) & (pivGoodData.UTM_X > 422000) \
                   & (pivGoodData.UTM_Y > 4608000) & (pivGoodData.UTM_Y < 4634000) & (pivGoodData.StationType=='Well')]
princpiv.drop_duplicates(subset = ['SampleId'],inplace=True)

ResOldPrinc = resultsNoND[(resultsNoND.SampleId.isin(princpiv.SampleId))]

GWStat = Stat[Stat.StationType.isin(['Well','Spring'])]
GWRes = results[results.StationId.isin(list(GWStat.StationId))]

Nitrate = GWRes[GWRes['ParAbb'].isin(['N','NO2','NO3','NH4'])]
NitrateStat = GWStat[GWStat.StationId.isin(list(Nitrate.StationId))]

"""
Explanation: Subset Data
End of explanation
"""


ParrAbbSummary = ResOldPrinc.groupby('ParAbb')['ResValue'].agg({'min':np.min, 'mean':np.mean,
                                                                'qrt5':(lambda x: np.percentile(x,q=5)),
                                                                'qrt95':(lambda x: np.percentile(x,q=95)),
                                                                'range':(lambda x: np.max(x)-np.min(x)),   
                                                                'lqrt':(lambda x: np.percentile(x,q=25)),
                                                                   'median':np.median, 
                                                                   'uqrt':(lambda x: np.percentile(x,q=75)), 
                                                                   'max':np.max, 'std':np.std, 
                                                                   'cnt':(lambda x: np.count_nonzero(~np.isnan(x)))}).reset_index()

ParrAbbSummary

manyPars = list(ParrAbbSummary[ParrAbbSummary['cnt'] >= 30]['ParAbb'])

ResOldPrinc = ResOldPrinc[ResOldPrinc['ParAbb'].isin(manyPars)] 
summaryStats = ParrAbbSummary[ParrAbbSummary['ParAbb'].isin(manyPars)]

summaryStats

from pylab import rcParams
rcParams['figure.figsize'] = 15, 10


parLabCounts = ParrAbbSummary.reset_index()
parLabCounts = parLabCounts.set_index(['ParAbb'])
plt.figure()
boxres= ResOldPrinc[ResOldPrinc['ParAbb'].isin(['pH_lab','pH_field'])] 
boxres.boxplot(column='ResValue', by='ParAbb',vert=False)

plt.title('Boxplot of Principal Aquifer pH')
plt.yticks([1,2],['Field pH (n = %s)'%(parLabCounts.loc['pH_field','cnt']),'Lab pH (n = %s)'%(parLabCounts.loc['pH_lab','cnt'])])
plt.xlim(6,9)
plt.xticks(np.arange(6,9.25,0.25))
plt.xlabel('pH')
plt.savefig(rootname+"pHBoxplot.svg")
plt.savefig(rootname+"pHBoxplot.pdf")

plt.figure()
boxres= ResOldPrinc[ResOldPrinc['ParAbb'].isin(['Temp'])] 
boxres.boxplot(column='ResValue', by='ParAbb',vert=False)

plt.title('Boxplot of Principal Aquifer Temperature')
plt.yticks([1],['Temperature (deg. C) (n = %s)'%(parLabCounts.loc['Temp','cnt'])])
plt.xticks(np.arange(5,30,1))
plt.xlabel('Temp. (deg. C)')
plt.savefig(rootname+"pHBoxplot.pdf")

import seaborn as sns
import matplotlib.pyplot as plt
sns.set(style="whitegrid")
rcParams['figure.figsize'] = 15, 20


parLabCounts = ParrAbbSummary.reset_index()
parLabCounts = parLabCounts.set_index(['ParAbb'])

parlist = ['Mg','Ca','Na','Cl','SO4','HCO3','Si','K','NO3','TDS','N']
boxres = ResOldPrinc[ResOldPrinc['ParAbb'].isin(parlist)] 
plt.figure()
sns.violinplot(x="ResValue", y='ParAbb', data=boxres, palette="Set3", scale='width', cut=0)
plt.xlabel('mg/L')
plt.xlim(0,1200)
plt.ylabel('Chemical Constituent')
plt.savefig(rootname+'violinMajor.pdf')

parLabCounts = ParrAbbSummary.reset_index()
parLabCounts = parLabCounts.set_index(['ParAbb'])


def parboxplot(parlist):
    plt.figure()
    boxres= ResOldPrinc[ResOldPrinc['ParAbb'].isin(parlist)] 
    boxres.boxplot(column='ResValue', by='ParAbb',vert=False)
    #labs = [str(parlist[i]) + " (n= %s)"%(parLabCounts.loc[parlist[i],'cnt']) for i in range(len(parlist))]
    #tickloc = [b+1 for b in range(len(parlist))]
    #plt.yticks(tickloc,labs)
    

parlist = ['pH_lab','pH_field']
parboxplot(parlist)
plt.xlabel('pH')
plt.savefig(rootname+'pHBoxplot')

parlist = ['Mg','Ca','Na','Cl','SO4','HCO3','Si','K','NO3','TDS','N']
parboxplot(parlist)
plt.title('Major Ions')
plt.xlabel('mg/L')
plt.grid(which='both',axis='both')
plt.xscale('log')
plt.xlim(0.1,1000)

plt.savefig(rootname+'MajorIonsBoxplot.pdf')

#plt.xlim(0.00001,1000)
#plt.xscale('log')

"""
Explanation: Summarize & Plot Data
End of explanation
"""


pipergrps.to_csv(rootname+'avgpiper.csv',index_label='StationId')
pivdata.to_csv(rootname+'pivotdata.csv',index_label='OBJECTID')
princpiv.to_csv(rootname+'PrincAquiferData.csv',index_label='OBJECTID')
pivgrp.to_csv(rootname+'pivgrps.csv',index_label='StationId')

Nitrate.to_csv(rootname+'NitrateResults.csv')
NitrateStat.to_csv(rootname+'NitrateStations.csv')

summaryStats.to_csv(rootname+'PrincAquifStats.csv')

summaryStats.to_clipboard()

GWStat.to_csv(rootname+'GWStations.csv',index_label='ObjectID')
GWRes.to_csv(rootname+'GWResults.csv',index_label='ObjectID')

writer = pd.ExcelWriter(rootname + "combined_out.xlsx", engine="xlsxwriter")
Stat.to_excel(writer, "stations", index=False)
results.to_excel(writer, "results", index=False)
GWStat.to_excel(writer, 'GWStations',index=False)
GWRes.to_excel(writer, 'GWResults',index=False)
pipergrps.to_excel(writer,'avgpiper')
pivdata.to_excel(writer,'pivotdata')
summaryStats.to_excel(writer,'princaquifstats')
writer.save()

"""
Explanation: Export data
End of explanation
"""


