## This code reclassifies the National Land Cover Dataset (NLCD) (30 m resolution) in the area of overlap
## with the 100% consensus MaxEnt model within the study area for a species with values for the Land Cover Utilization Index
## calculated by the R program "NLCDPercentCoveriNat15mMaxEntFSELoop.R"

## NOTE: Use ArcMap project template LandCoverSelectivityAssessmentMap_AoenNoLegends.mxd

## NOTE: The geoprocessing environments will be set below in python script 

# Make arcmap python commands and spatial analyst (sa) operations available
import arcpy
from arcpy.sa import *

import os

# Overwrite pre-existing files
arcpy.env.overwriteOutput = True

################################################  
# set working directory
arcpy.env.workspace = "K:/Monarch"

# Set input and output directories
Species = "A_cap"
Species2 = "Acap"
InDirect = "K:/Monarch/" 
OutDirect = "K:/Monarch/"
MxtCoreName = Species2 + "100con1k" # Name of MaxEnt Core model (100% consensus FSE)
MxtCoreName2 = Species2 + "mxtcor"
LandCoverRasterName = "nlcd2011scsta"
StudyAreaPolyName = "SCStatesEcoregionsDissAlbers.shp"
LandAreaPolyName = "USGSAdminLandAlbers.shp"

import csv
with open(InDirect + Species + "_ReClassifyNLCDtoRelativizedElectivityIndexPreferred.csv") as f:
    reader = csv.reader(f)
    next(reader) # skip header
    ReClassifyDataPreferred = [r for r in reader]

import csv
with open(InDirect + Species + "_ReClassifyNLCDtoRelativizedElectivityIndexNeutral.csv") as f:
    reader = csv.reader(f)
    next(reader) # skip header
    ReClassifyDataNeutral = [r for r in reader]

import csv
with open(InDirect + Species + "_ReClassifyNLCDtoRelativizedElectivityIndexAvoided.csv") as f:
    reader = csv.reader(f)
    next(reader) # skip header
    ReClassifyDataAvoided = [r for r in reader]


# Project MaxentRasterCore grid to Albers
arcpy.ProjectRaster_management(in_raster=InDirect + MxtCoreName, out_raster=InDirect + MxtCoreName2 + "A",
                               out_coor_system="PROJCS['North_America_Albers_Equal_Area_Conic',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-96.0],PARAMETER['Standard_Parallel_1',20.0],PARAMETER['Standard_Parallel_2',60.0],PARAMETER['Latitude_Of_Origin',40.0],UNIT['Meter',1.0]]", resampling_type="NEAREST", cell_size="906.71730152707 906.71730152707", geographic_transform="", Registration_Point="", in_coor_system="PROJCS['NAD_1983_Albers',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-96.0],PARAMETER['Standard_Parallel_1',29.5],PARAMETER['Standard_Parallel_2',45.5],PARAMETER['Latitude_Of_Origin',23.0],UNIT['Meter',1.0]]", vertical="NO_VERTICAL")

# In ArcMap project, set Geoprocessing Environments to have snap raster and cell size equivalent to that of above Albers Maxent core raster
arcpy.env.snapRaster = Raster(InDirect + MxtCoreName2 + "A")
arcpy.env.cellSize = Raster(InDirect + MxtCoreName2 + "A")
# Set processing extent to land area
arcpy.env.extent = InDirect + LandAreaPolyName

# Set to null values < 1 in above raster
arcpy.gp.SetNull_sa(InDirect + MxtCoreName2 + "A", "1", InDirect + MxtCoreName2 + "A1", "VALUE < 1")

# Clip above raster to study area
arcpy.Clip_management(in_raster=InDirect + MxtCoreName2 + "A1",
                      rectangle="-645169.704274798 -1645583.74593679 655251.099961824 -331982.678015766",
                      out_raster=InDirect + MxtCoreName2 + "As", in_template_dataset=InDirect + StudyAreaPolyName,
                      nodata_value="0", clipping_geometry="ClippingGeometry", maintain_clipping_extent="NO_MAINTAIN_EXTENT")

# Read in MaxEnt core raster clipped to study area
MxtCoreStudyArea = Raster(InDirect + MxtCoreName2 + "As")

# Create zero version of above raster
MxtCoreStudyZero = ((MxtCoreStudyArea + 1)/(MxtCoreStudyArea + 1))-1
MxtCoreStudyZero.save(InDirect + MxtCoreName2 + "As0")

# In ArcMap project, set Geoprocessing Environments to have snap raster and cell size equivalent to that of the NLCD raster
arcpy.env.snapRaster = Raster(InDirect + LandCoverRasterName)
arcpy.env.cellSize = Raster(InDirect + LandCoverRasterName)
# Set processing extent to study area
arcpy.env.extent = InDirect + StudyAreaPolyName

# Restrict NLCD raster by adding above zero version of Maxent Core Raster clipped to study area to make
# base for Relative Electivity Index (RELELE)
# MxtCoreStudyZero=Raster(InDirect + MxtCoreName2 + "As0")
MxtNLCDRaster = Plus(LandCoverRasterName, MxtCoreStudyZero)
MxtNLCDRaster.save(OutDirect+Species2+"RELELE")

# Reclassify above raster for Preferred Land cover using ReClassifyDataPreferred csv file values (integer relativized electivity x 10,0000) for land cover classes to create preferred raster
# (omitted classes, such as water, will be assigned to No Data)
LandCoverPreferredRaster1 = Reclassify(OutDirect+Species2+"RELELE", "VALUE", RemapValue(ReClassifyDataPreferred), "NODATA")

# Convert above raster back to floating point with decimal values of Relativized Electivity
LandCoverPreferredRaster2 = Float(LandCoverPreferredRaster1)/10000
LandCoverPreferredRaster2.save(OutDirect+Species2+"RELPREF")

# Reclassify above raster for Avoided Land cover using ReClassifyDataAvoided csv file values (integer relativized electivity x 10,0000) for land cover classes to create preferred raster
# (omitted classes, such as water, will be assigned to No Data)
LandCoverAvoidedRaster1 = Reclassify(OutDirect+Species2+"RELELE", "VALUE", RemapValue(ReClassifyDataAvoided), "NODATA")

# Convert above raster back to floating point with decimal values of Relativized Electivity
LandCoverAvoidedRaster2 = Float(LandCoverAvoidedRaster1)/10000
LandCoverAvoidedRaster2.save(OutDirect+Species2+"RELAVOID")


# Reclassify above raster for Neutral Land cover using ReClassifyDataNeutral csv file values (integer relativized electivity x 10,0000) for land cover classes to create preferred raster
# (omitted classes, such as water, will be assigned to No Data)
LandCoverNeutralRaster1 = Reclassify(OutDirect+Species2+"RELELE", "VALUE", RemapValue(ReClassifyDataNeutral), "NODATA")

# Convert above raster back to floating point with decimal values of Relativized Electivity
LandCoverNeutralRaster2 = Float(LandCoverNeutralRaster1)/10000
LandCoverNeutralRaster2.save(OutDirect+Species2+"RELANEUT")

####################################################################
### Display the three output rasters in ArcMap using stretched color ramps,
### Refer to the example template of the ArcMap project
### "LandCoverRiskAssessmentElectivityMap_AoenLegends.mxd"
#####################################################################
