# SUMMARY:      rowcolmap.py
# USAGE:        Part of python version createstreamnetwork. 
#               Original Description: It is a simple script used to produce a
#               polygon representation of a grid. Each polygon in outcover
#               represents a single cell, and is assigned row and column 
#               numbers.
# ORG:          Pacific Northwest National Laboratory
# E-MAIL:       zhuoran.duan@pnnl.gov
# ORIG-DATE:    Apr-2017
# DESCRIPTION:  Python version of original createstreamnetwork aml
# DESCRIP-END.
# COMMENTS:     This python script is created based on original 
#		AML scripts roadaspect.aml as part of DHSVM.
#
# Last Change: 2017-08-10 by zduan


import arcpy
from arcpy import env
from arcpy.sa import *
import arcgisscripting
import numpy as np

def rowcolmapfun(elev, rcpoly, maskExt):
	demraster = arcpy.Raster(elev)
	env.mask = maskExt

	# Get raster properties
	rows = int(arcpy.GetRasterProperties_management(demraster, "ROWCOUNT")[0])
	cols = int(arcpy.GetRasterProperties_management(demraster, 'COLUMNCOUNT')[0])
	csizex = float(arcpy.GetRasterProperties_management(demraster, 'CELLSIZEX')[0])
	csizey = float(arcpy.GetRasterProperties_management(demraster, 'CELLSIZEY')[0])

	# Create a dummy raster
	arr = np.arange(rows * cols).reshape(rows, cols)
	lowerLeft = arcpy.Point(demraster.extent.XMin, demraster.extent.YMin)
	dummy = arcpy.NumPyArrayToRaster(arr, lowerLeft, csizex, csizey)
	
	# Convert to feature
	if arcpy.Exists(rcpoly):
		arcpy.Delete_management(rcpoly)

	arcpy.RasterToPolygon_conversion(dummy, rcpoly, "NO_SIMPLIFY", "VALUE")
	arcpy.DefineProjection_management(rcpoly, demraster.extent.spatialReference)
	arcpy.AddField_management(rcpoly, "ROW", "LONG")
	arcpy.AddField_management(rcpoly, "COL", "LONG")

	# Update columns and rows
	targfields = ['GRIDCODE', 'COL', 'ROW']
	with arcpy.da.UpdateCursor(rcpoly, targfields) as recs:
		for rec in recs:
			rec[1] = rec[0] % cols
			rec[2] = rec[0] / cols
			recs.updateRow(rec)

	print('Join table finished')

	arcpy.Delete_management(dummy)