# ---------------------------------------------------------------------------
# extract_basin_prism_values.py
# Created on: 2014-07-22 18:30:25.00000 (generated by ArcGIS/ModelBuilder)
# Description: extract prism data from raster using basin shapefiles as mask
# and output a csv file for each basin
# UPDATED (8/4/2014): calculate basin mean temperature or precipitation
# Created by Ryan Spies (ryan.spies@amec.com)
# ---------------------------------------------------------------------------

################### User Input #####################
RFC = 'WGRFC_FY2015'
variable = 'ppt' # use temperature: 'tmean' or precipitation: 'ppt'
basins_folder = r'P:\NWS\GIS\WGRFC_FY2015\Shapefiles\calb_basins'

# if you only want to run specific basins -> list them below
# otherwise set it equal to empty list (basins_overwrite = [])
basins_overwrite = [] 

output_dir = 'P:\\NWS\\GIS\\' + RFC + '\\PRISM\\Monthly_' + variable +'\\'  
################# End User Input ######################
print '#### Make sure these directory exists:\nC:\\NWS\\python\\temp_output\\ \n' + output_dir
raw_input("Press enter to begin processing...")

# Years to create PRISM Raster table (CONUS)
#years = [1994,1995,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010]
years = [2001,2002,2003,2004,2005,2006,2007,2008,2009,2010]

# Import arcpy module
import arcpy
import os
import csv
arcpy.env.overwriteOutput = True

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "P:\\NWS\\GIS\\Models\\10_0_tools\\Model_Output.gdb" # temporary file storage directory
arcpy.env.parallelProcessingFactor = "50"
print 'ok so far...'

#################################################################################
# find all basins in RFC or only run the specified basin list
basin_files = os.listdir(basins_folder) # list all basin shapefiles in the above specified directory
if len(basins_overwrite) != 0:
    basin_files = basins_overwrite # use the basins_overright variable to run only specified basins instead of all RFC basins
basins = []
check_dir = os.listdir(output_dir) # list all folders in output_dir
for each in basin_files:
    if each[:5] not in basins:
        basins.append(each[:5])

print 'Identified ' + str(len(basins)) + ' basins in ' + RFC + ' input directory...'        
if variable == 'ppt':
    label = ' (mm) '
if variable == 'tmean':
    label = ' (F) '
# loop through basins

for basin in basins:
    ## Script arguments
    Basin_Boundary = basins_folder + '\\' + basin + '.shp'
    prism_csv = open(output_dir + basin + '_prism' + '.csv', 'wb')
    csvFile = csv.writer(prism_csv) #output csv
    csvFile.writerow(['year','month',variable + label])
            
    print basin 
    for year in years:
        month = 1
        while month <= 12:
            print 'month:' + str(year) + ' ' + str(month)
            if month < 10:
                mncount = '0' + str(month)
            else:
                mncount = str(month)
            PRISM_Dataset = 'Q:\\GISLibrary\\PRISM\\Calendar_Year\\' + variable + '\\'+str(year)+'\\' + '\\PRISM_'+variable+'_stable_4kmM2_'+str(year)+mncount+'_asc.asc'

            ## Local variables: Changed to output to C: drive (errors out when on the P:)
            #Basin_Raster = 'P:\\NWS\\GIS\\Models\\10_0_tools\\Model_Output\\' + basin
            #Basin_Points = 'P:\\NWS\\GIS\\Models\\10_0_tools\\Model_Output\\' + basin + '_points'
            #Stats_Table = 'P:\\NWS\\GIS\\Models\\10_0_tools\\Model_Output\\prism_stats.dbf'
            Basin_Raster = 'C:\\NWS\\python\\temp_output\\' + basin
            Basin_Points = 'C:\\NWS\\python\\temp_output\\' + basin + '_points'
            Stats_Table = 'C:\\NWS\\python\\temp_output\\prism_stats.dbf'
        
            ## Process: Extract by Mask
            print 'Extracting by mask...'
            arcpy.gp.ExtractByMask_sa(PRISM_Dataset, Basin_Boundary, Basin_Raster)
        
            ## Process: Raster to Point
            print 'Raster to point...'
            arcpy.RasterToPoint_conversion(Basin_Raster, Basin_Points, "VALUE")
            
            # Process: Summary Statistics
            print 'Calculating Summary statistics...'
            arcpy.Statistics_analysis(Basin_Points + '.shp', Stats_Table, "GRID_CODE MEAN", "")
            
            # Process: output csv file
            print 'Creating '+ basin + '_prism.csv file...'
            rows = arcpy.SearchCursor(Stats_Table)
            
            fieldnames = [f.name for f in arcpy.ListFields(Stats_Table)]
        
            allRows = []
            for row in rows:
                rowlist = []
                for field in fieldnames:
                    rowlist.append(row.getValue(field))
                    print rowlist
                if variable != 'ppt':
                    final_var=((float(rowlist[2])*1.8)+32)
                else:
                    final_var=float(rowlist[2])
                print variable + ': ' + str(final_var)
            csvFile.writerow([year,month,final_var])
            month += 1

    prism_csv.close()
    
print 'Completed grid extraction!'
