# -*- coding: cp1252 -*-
## This code converts spatially thinned high resolution (app. 30 m) presence points shapefile to North America Albers Equal Area Conic
## projection and add nearness to roads. 

## NOTE: If available, can use ArcMap project template LandCoverRiskAssessment_PointProcessing.mxd

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


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/"
OutDirect2 = "K:/Monarch"
NLCDYear = "2016"
SpatFiltBuff = "1" # spatial filter in km
NAAlbersPoly = "AdminGlobalAlbers.shp" # NA Albers Equal Area Conic projected shapefile
SCStatesAreaPolyName = "SCStatesEcoregionsDissGeo.shp"
RoadPolylineName = "OKARTXLA_OSM_Roads_alb"
NLCDRasterName = "nlcd2016scst" 
MxtCoreName = Species2 + "100con1k" # Name of MaxEnt Core model (100% consensus FSE) with zero values as null
MxtCoreSCStatesPolyNameAlbers = Species2 + "MxtCoreSCPolyAlbers.shp"
MxtCoreSCStatesPolyNameWGS84 = Species2 + "MxtCoreSCPolyWGS84.shp"
PresencePointsUnClippedWGS84 = Species + "_HiRes_UnClipped_Unthinned.shp"
TargetGroupPointsUnClippedWGS84 = Species + "_TargetGroupHiRes_UnClipped_Unthinned.shp"

#####################################################################################
### PART 1: Use intersection of core MaxEnt model and SCStates to define study area and tabulate area of NLCD classses for 
### study area 
#########################################################################################
#
### Read in MaxEnt Core habitat (100% Consensus for Subset Ensemble (FSE))
MxtCore = Raster(InDirect + MxtCoreName)

# In ArcMap project, set Geoprocessing Environments to have snap raster and cell size equivalent to that of MaxEnt FSE
arcpy.env.outputCoordinateSystem = InDirect + NAAlbersPoly
arcpy.env.snapRaster = MxtCore
arcpy.env.cellSize = MxtCore

# Set Null zero values in Maxent Core
#arcpy.gp.SetNull_sa(InDirect + MxtCoreName, InDirect + MxtCoreName, InDirect + MxtCoreName + "1", "VALUE < 1")

# Clip above raster to SCStatesArea
arcpy.Clip_management(in_raster=InDirect + MxtCoreName,  out_raster=InDirect + MxtCoreName + "C",
                      in_template_dataset=InDirect + SCStatesAreaPolyName,
                      nodata_value="0", clipping_geometry="ClippingGeometry", maintain_clipping_extent="NO_MAINTAIN_EXTENT")

# Convert above raster to polygon
arcpy.RasterToPolygon_conversion(in_raster=InDirect + MxtCoreName + "C", out_polygon_features=InDirect+MxtCoreSCStatesPolyNameAlbers, simplify="NO_SIMPLIFY", raster_field="VALUE",
                                 create_multipart_features="SINGLE_OUTER_PART", max_vertices_per_feature="")

# Project above polygon shapefile to wgs84
arcpy.Project_management(in_dataset=InDirect+MxtCoreSCStatesPolyNameAlbers, out_dataset=InDirect+MxtCoreSCStatesPolyNameWGS84,
                         out_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", transform_method="WGS_1984_(ITRF00)_To_NAD_1983", in_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]]", preserve_shape="NO_PRESERVE_SHAPE", max_deviation="", vertical="NO_VERTICAL")

# Project clipped Maxent core raster to North American Albers Equal Area Conic
#arcpy.ProjectRaster_management(in_raster=InDirect + MxtCoreName + "C", out_raster=InDirect + MxtCoreName + "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="977.887740603613 977.887740603613", geographic_transform="WGS_1984_(ITRF00)_To_NAD_1983", Registration_Point="", in_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")

# Use Tabulate Area tool to obtain areas of all land cover types for the above MaxEnt raster
arcpy.gp.TabulateArea_sa(InDirect + MxtCoreName + "C", "VALUE", InDirect + NLCDRasterName,
                         "VALUE", InDirect + MxtCoreName + "_" + NLCDYear + "nlcdstats.dbf", InDirect + MxtCoreName + "C")


############################################################################################
### PART 2: Clip to MaxEnt Core Habitat within Study Area and spatially thin high resolution
### iNaturalist presence and target group point data processed in R script
############################################################################################
#
## First clip presence shapefile to Maxent Core Habitat within Study Area
arcpy.Clip_analysis(in_features=InDirect + PresencePointsUnClippedWGS84, clip_features=InDirect+MxtCoreSCStatesPolyNameWGS84,
                    out_feature_class=InDirect + Species + "_HiRes_Unthinned.shp", cluster_tolerance="")

## Then spatially thin to 1 km
## Create buffer around points using SpatFiltBuff
arcpy.Buffer_analysis(in_features=InDirect + Species + "_HiRes_Unthinned.shp",
                      out_feature_class=InDirect + Species + "_HiRes_Unthinned_" + SpatFiltBuff + "kmBuff.shp",
                      buffer_distance_or_field=SpatFiltBuff, line_side="FULL", line_end_type="ROUND",
                      dissolve_option="NONE", dissolve_field="", method="PLANAR")

## Delete identical in buffer using SpatFiltBuff
arcpy.DeleteIdentical_management(in_dataset=InDirect + Species + "_HiRes_Unthinned_" + SpatFiltBuff + "kmBuff.shp", fields="Shape",
                                 xy_tolerance=SpatFiltBuff, z_tolerance="0")

## Convert buffer file to points to get spatially thinned data
arcpy.FeatureToPoint_management(in_features=InDirect +  Species + "_HiRes_Unthinned_" + SpatFiltBuff + "kmBuff.shp",
                                out_feature_class=InDirect + Species + "_HiRes_Thinned" + SpatFiltBuff + "km.shp", point_location="CENTROID")

### Clip Target Group shapefile to above convex hull polygon of thinned presence points

# First Create convex hull polygon for presence points
# Create convex hull polygon around Albers presence points in study area
arcpy.MinimumBoundingGeometry_management(in_features=OutDirect + Species + "_HiRes_Thinned1km.shp",
                                         out_feature_class=OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly.shp", geometry_type="CONVEX_HULL", group_option="ALL", group_field="", mbg_fields_option="NO_MBG_FIELDS")

# Create 20 km buffer around above convex hull polygon
arcpy.Buffer_analysis(in_features=OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly.shp",
                      out_feature_class=OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly20kBuffer1.shp", buffer_distance_or_field="20 Kilometers", line_side="FULL", line_end_type="ROUND", dissolve_option="NONE", dissolve_field="", method="GEODESIC")

# Clip above buffer to study area
arcpy.Clip_analysis(in_features=OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly20kBuffer1.shp",
                    clip_features=InDirect + SCStatesAreaPolyName,
                    out_feature_class=OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly20kBuffer.shp", cluster_tolerance="")

## Now clip target group shapefile to convex hull polygon
arcpy.Clip_analysis(in_features=InDirect + TargetGroupPointsUnClippedWGS84, clip_features=InDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly20kBuffer.shp",
                    out_feature_class=InDirect + Species + "_TargetGroupHiRes_Unthinned.shp", cluster_tolerance="")

## Then spatially thin to 1 km
## Create buffer around points using SpatFiltBuff
arcpy.Buffer_analysis(in_features=InDirect + Species + "_TargetGroupHiRes_Unthinned.shp",
                      out_feature_class=InDirect + Species + "_TargetGroupHiRes_Unthinned_" + SpatFiltBuff + "kmBuff.shp",
                      buffer_distance_or_field=SpatFiltBuff, line_side="FULL", line_end_type="ROUND", dissolve_option="NONE",
                      dissolve_field="", method="PLANAR")

## Delete identical in buffer using SpatFiltBuff
arcpy.DeleteIdentical_management(in_dataset=InDirect + Species + "_TargetGroupHiRes_Unthinned_" + SpatFiltBuff + "kmBuff.shp",
                                 fields="Shape", xy_tolerance=SpatFiltBuff, z_tolerance="0")

## Convert buffer file to points to get spatially thinned data
arcpy.FeatureToPoint_management(in_features=InDirect +  Species + "_TargetGroupHiRes_Unthinned_" + SpatFiltBuff + "kmBuff.shp",
                                out_feature_class=InDirect + Species + "_TargetGroupHiRes_Thinned" + SpatFiltBuff + "km.shp", point_location="CENTROID")

## Sample 10,000 out of all above points
arcpy.CreateRandomPoints_management(out_path=OutDirect2, out_name=Species + "_TargetGroupHiRes_Thinned1km_10000a",
                                    constraining_feature_class=Species + "_TargetGroupHiRes_Thinned1km",
                                    number_of_points_or_field="10000",
                                    minimum_allowed_distance="0 Meters", create_multipoint_output="POINT", multipoint_size="0")

## Intersect new points with original data
arcpy.Intersect_analysis(in_features=[OutDirect + Species + "_TargetGroupHiRes_Thinned1km.shp", OutDirect + Species + "_TargetGroupHiRes_Thinned1km_10000a.shp"],
                         out_feature_class=OutDirect + Species + "_TargetGroupHiRes_Thinned1km_10000.shp",
                         join_attributes="ALL", cluster_tolerance="-1 Unknown", output_type="INPUT")

############################################################################################
### PART 3: Add distances to roads to thinned high resolution iNaturalist presence and target group point data processed in R script
############################################################################################

## Process presence points

# Convert clipped presence point data to Albers

# Use Near Tool to add distance to nearest road. Can take about 20 minutes.
arcpy.Near_analysis(in_features=OutDirect + Species + "_HiRes_Thinned1km.shp",
                    near_features=InDirect + RoadPolylineName + ".shp", search_radius="", location="NO_LOCATION", angle="NO_ANGLE", method="GEODESIC")

# Project above shapefile back to Geographic WGS84
arcpy.Project_management(in_dataset=OutDirect + Species + "_HiRes_Thinned1km.shp",
                         out_dataset=OutDirect + Species + "_HiRes_Thinned1km_rd.shp",
                         out_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],VERTCS['Unknown VCS',VDATUM['Unknown'],PARAMETER['Vertical_Shift',0.0],PARAMETER['Direction',1.0],UNIT['Meter',1.0]]", transform_method="WGS_1984_(ITRF00)_To_NAD_1983", in_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]]", preserve_shape="NO_PRESERVE_SHAPE", max_deviation="")

###

## Process Target Group

# Convert clipped presence point data to Albers

# Use Near Tool to add distance to nearest road. Can take about 20 minutes.
arcpy.Near_analysis(in_features=OutDirect + Species + "_TargetGroupHiRes_Thinned1km_10000.shp",
                    near_features=InDirect + RoadPolylineName + ".shp", search_radius="", location="NO_LOCATION", angle="NO_ANGLE", method="GEODESIC")

# Project above shapefile back to Geographic WGS84
arcpy.Project_management(in_dataset=OutDirect + Species + "_TargetGroupHiRes_Thinned1km_10000.shp",
                         out_dataset=OutDirect + Species + "_TargetGroupHiRes_Thinned1km_10000rd.shp",
                         out_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],VERTCS['Unknown VCS',VDATUM['Unknown'],PARAMETER['Vertical_Shift',0.0],PARAMETER['Direction',1.0],UNIT['Meter',1.0]]", transform_method="WGS_1984_(ITRF00)_To_NAD_1983", in_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]]", preserve_shape="NO_PRESERVE_SHAPE", max_deviation="")


### Next, create 10,000 random points for buffered convex hull polygon around high resolution
### presence points within study area and add distances to roads and later addin 2015 population density per 1 km in R


# Create 10,000 random points within above polygon with minimum 1 km distance
arcpy.CreateRandomPoints_management(out_path=OutDirect, out_name=Species + "_HiRes_ConvexHull_10000points",
                                    constraining_feature_class=OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly20kBuffer.shp", constraining_extent="0 0 250 250", number_of_points_or_field="10000", minimum_allowed_distance="1 Kilometers", create_multipoint_output="POINT", multipoint_size="0")

# Use Near Tool to add distance to nearest road. Can take about 40 minutes
# First, set processing extent to buffered convex hull area
arcpy.env.extent = OutDirect + Species + "_HiRes_Thinned1km_ConvexHullPoly20kBuffer.shp"
arcpy.Near_analysis(in_features=OutDirect + Species + "_HiRes_ConvexHull_10000points.shp", near_features=InDirect + RoadPolylineName + ".shp", search_radius="", location="NO_LOCATION", angle="NO_ANGLE", method="GEODESIC")

# Project above shapefile to Geographic WGS84
arcpy.Project_management(in_dataset=OutDirect + Species + "_HiRes_ConvexHull_10000points.shp",
                         out_dataset=OutDirect + Species + "_HiRes_ConvexHull_10000pointsg.shp", out_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],VERTCS['Unknown VCS',VDATUM['Unknown'],PARAMETER['Vertical_Shift',0.0],PARAMETER['Direction',1.0],UNIT['Meter',1.0]]", transform_method="WGS_1984_(ITRF00)_To_NAD_1983", in_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]]", preserve_shape="NO_PRESERVE_SHAPE", max_deviation="")

############################################################

#############################################################################################
### Return to R script "LandCoverRiskAssessmentRelativeChange_NLCDPercentCoverHiResVsRandom.r" to
### weight presence points accoring to population and distance to road and calculate
### land cover utilization indices for each MaxEnt model
#############################################################################################
