# ---------------------------------------------------------------------------
# extract_basin_DEM_statistics.py
# Created on: 2014-07-22 18:31:56.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: extract_basin_DEM_statistics <Basin> 
# Description: extract elevation and slope characteristics for a basin shapefile
# using the user define DEM raster
# Created and Modified by Ryan Spies
# ---------------------------------------------------------------------------
print 'Importing modules...'
# Import arcpy module
import arcpy
import os
import csv
arcpy.env.overwriteOutput = True
os.chdir("../..")
maindir = os.getcwd()

################### User Input #####################
RFC = 'MBRFC_FY2015'
basins_folder = maindir + '\\' + RFC[:5] + os.sep + RFC + '\\Shapefiles\\calb_basins\\'
dem_file = maindir + '\\GTOPO_1K_Hydro\\na_dem_null' # 1km resolution
#dem_file = r'Q:\GISLibrary\NHDPlus\v1\WG\NHDPlus12\Elev_Unit_d\elev_cm'
#dem_file = r'P:\NWS\GIS\APRFC\Processing\ASTER\kusk_mosiac.tif'

output_dir = maindir + '\\'+ RFC[:5] + os.sep + RFC + '\\Elevation_Slope\\Stats_Out\\'
sr = arcpy.SpatialReference(4269) # define projection of basin shp -> 4269 = GCS_North_American_1983

# manually enter basins to analyze (example: basins = ['MFDN3'])
# otherwise leave blank and script will
# analyze all basins in the basins_folder specified above
basins_overwrite = [] 
################# End User Input ######################

if not os.path.exists('C:\\NWS\\python\\temp_output\\'):
    print "Missing directory: 'C:\\NWS\\python\\temp_output\\' -> please create"
    raw_input("Press enter to continue processing...")
if not os.path.exists(output_dir):
    print "Missing directory: " + output_dir + " -> please create"
    raw_input("Press enter to continue processing...")

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

# Local variables (temporary files):
Basin_DEM = 'C:\\NWS\\python\\temp_output\\basin_DEM'
Slope_Raster = 'C:\\NWS\\python\\temp_output\\basin_slope'
Slope_Points = 'C:\\NWS\\python\\temp_output\\slope_points'
Slope_Stats__dbf = 'C:\\NWS\\python\\temp_output\\slope_stats.dbf'
Elevation_Points = 'C:\\NWS\\python\\temp_output\\elevation_pts'
Elevation_Stats__dbf = 'C:\\NWS\\python\\temp_output\\elevation_stats'

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "C:\\NWS\\python\\Model_Output.gdb" # temporary file storage directory
#arcpy.env.parallelProcessingFactor = "50"
print 'ok so far...'

# find all basins in RFC task or only run the specified basin overwrite 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.split('.')[0] not in basins:
        basins.append(each.split('.')[0])
            
for basin in basins:    
    print basin
    Basin_Boundary = basins_folder + '\\' + basin + '.shp'
    check_project = basins_folder + '\\' + basin + '.prj'

    # Process: Define Projection
    print 'Define Projection...'
    if not os.path.exists(check_project):
        arcpy.DefineProjection_management(Basin_Boundary, sr)
    
    # Process: Extract by Mask
    print 'Extract by mask...'
    arcpy.gp.ExtractByMask_sa(dem_file, Basin_Boundary, Basin_DEM)
    
    # Process: Slope
    print 'Calculating basin slope raster...'
    arcpy.gp.Slope_sa(Basin_DEM, Slope_Raster, "PERCENT_RISE", "0.01")
    
    # Process: Raster to Point (2)
    print 'Coverting raster to point...'
    arcpy.RasterToPoint_conversion(Slope_Raster, Slope_Points, "Value")
    
    # Process: Summary Statistics (2)
    print 'Calculating mean slope...'
    arcpy.Statistics_analysis(Slope_Points + '.shp', Slope_Stats__dbf, "grid_code MEAN", "")

    # Process: output csv file
    print 'Creating '+ basin + '_mean_slope_percent.csv file...'
    rows = arcpy.SearchCursor(Slope_Stats__dbf)
    slope_csv = open(output_dir + basin + '_mean_slope_percent.csv', 'wb')
    csvFile = csv.writer(slope_csv) #output csv
    fieldnames = [f.name for f in arcpy.ListFields(Slope_Stats__dbf)]

    allRows = []
    for row in rows:
        rowlist = []
        for field in fieldnames:
            rowlist.append(row.getValue(field))
        allRows.append(rowlist)

    csvFile.writerow(fieldnames)
    for row in allRows:
        csvFile.writerow(row)
    row = None
    rows = None
    slope_csv.close()
        
    # Process: Raster to Point
    print 'Converting elevation raster to point...'
    arcpy.RasterToPoint_conversion(Basin_DEM, Elevation_Points, "VALUE")
    
    # Process: Summary Statistics
    print 'Calculating elevation statistics...'
    arcpy.Statistics_analysis(Elevation_Points + '.shp', Elevation_Stats__dbf, "grid_code MIN;grid_code MEAN;grid_code MAX", "")
    
    # Process: output csv file
    print 'Creating '+ basin + '_elevation_stats_cm.csv file...'
    rows = arcpy.SearchCursor(Elevation_Stats__dbf)
    elevation_csv = open(output_dir + basin + '_elevation_stats_cm.csv', 'wb')
    csvFile = csv.writer(elevation_csv) #output csv
    fieldnames = [f.name for f in arcpy.ListFields(Elevation_Stats__dbf)]

    allRows = []
    for row in rows:
        rowlist = []
        for field in fieldnames:
            rowlist.append(row.getValue(field))
        allRows.append(rowlist)

    csvFile.writerow(fieldnames)
    for row in allRows:
        csvFile.writerow(row)
    row = None
    rows = None
    elevation_csv.close()
    
print 'Completed!!' 
