<jupyter_start><jupyter_text>Sweden  Post-Harvest stream buffers<jupyter_code>import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.parallelProcessingFactor = "100%"<jupyter_output><empty_output><jupyter_text>Select harvest sites that were harvested within 5 years before being scanned with LiDAR. Some areas were scanned multiple times. The shapefile 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/laserSquare_scandates.shp' contains a column named diff_years which indicates how many different years each tile was scanned with LiDAR. If the value is over 0 it means that it was scanned multiple years. Areas that were scanned multiple years were excluded to avoid confusion.<jupyter_code>swedish_lidar_original = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/laserSquare_scandates.shp'
scanned_one_year = arcpy.SelectLayerByAttribute_management(swedish_lidar_original, 'NEW_SELECTION', '"diff_years" = 0')
arcpy.CopyFeatures_management(scanned_one_year , 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp')

sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

## 2009
CC_5_years_before_2009 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2004 and "HARVYR" < 2009')
arcpy.CopyFeatures_management(CC_5_years_before_2009, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2004_2009.shp')

sweden_lidar_2009 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2009')
arcpy.CopyFeatures_management(sweden_lidar_2009, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2009.shp')

CC_2004_2009 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2004_2009.shp'
lidar_2009 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2009.shp'

harvested_5_years_before_lidar_2009 = arcpy.SelectLayerByLocation_management(CC_2004_2009, 'INTERSECT', lidar_2009,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2009, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2004_to_2009_lidar_2009.shp')


## 2010
CC_5_years_before_2010 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2005 and "HARVYR" < 2010')
arcpy.CopyFeatures_management(CC_5_years_before_2010, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2005_2010.shp')

sweden_lidar_2010 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2010')
arcpy.CopyFeatures_management(sweden_lidar_2010, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2010.shp')

CC_2005_2010 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2005_2010.shp'
lidar_2010 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2010.shp'

harvested_5_years_before_lidar_2010 = arcpy.SelectLayerByLocation_management(CC_2005_2010, 'INTERSECT', lidar_2010,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2010, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2005_to_2010_lidar_2010.shp')


## 2011
CC_5_years_before_2011 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2006 and "HARVYR" < 2011')
arcpy.CopyFeatures_management(CC_5_years_before_2011, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2006_2011.shp')

sweden_lidar_2011 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2011')
arcpy.CopyFeatures_management(sweden_lidar_2011, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2011.shp')

CC_2006_2011 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2006_2011.shp'
lidar_2011 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2011.shp'

harvested_5_years_before_lidar_2011 = arcpy.SelectLayerByLocation_management(CC_2006_2011, 'INTERSECT', lidar_2011,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2011, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2006_to_2011_lidar_2011.shp')

sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

## 2012
CC_5_years_before_2012 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2007 and "HARVYR" < 2012')
arcpy.CopyFeatures_management(CC_5_years_before_2012, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2007_2012.shp')

sweden_lidar_2012 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2012')
arcpy.CopyFeatures_management(sweden_lidar_2012, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2012.shp')

CC_2007_2012 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2007_2012.shp'
lidar_2012 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2012.shp'

harvested_5_years_before_lidar_2012 = arcpy.SelectLayerByLocation_management(CC_2007_2012, 'INTERSECT', lidar_2012,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2012, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2007_to_2012_lidar_2012.shp')

sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

## 2013
CC_5_years_before_2013 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2008 and "HARVYR" < 2013')
arcpy.CopyFeatures_management(CC_5_years_before_2013, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2008_2013.shp')

sweden_lidar_2013 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2013')
arcpy.CopyFeatures_management(sweden_lidar_2013, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2013.shp')

CC_2008_2013 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2008_2013.shp'
lidar_2013 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2013.shp'

harvested_5_years_before_lidar_2013 = arcpy.SelectLayerByLocation_management(CC_2008_2013, 'INTERSECT', lidar_2013,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2013, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2008_to_2013_lidar_2013.shp')

sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

## 2014
CC_5_years_before_2014 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2009 and "HARVYR" < 2014')
arcpy.CopyFeatures_management(CC_5_years_before_2014, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2009_2014.shp')

sweden_lidar_2014 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2014')
arcpy.CopyFeatures_management(sweden_lidar_2014, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2014.shp')

CC_2009_2014 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2009_2014.shp'
lidar_2014 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2014.shp'

harvested_5_years_before_lidar_2014 = arcpy.SelectLayerByLocation_management(CC_2009_2014, 'INTERSECT', lidar_2014,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2014, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2009_to_2014_lidar_2014.shp')

sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

## 2015
CC_5_years_before_2015 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2010 and "HARVYR" < 2015')
arcpy.CopyFeatures_management(CC_5_years_before_2015, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2010_2015.shp')

sweden_lidar_2015 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2015')
arcpy.CopyFeatures_management(sweden_lidar_2015, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2015.shp')

CC_2010_2015 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2010_2015.shp'
lidar_2015 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2015.shp'

harvested_5_years_before_lidar_2015 = arcpy.SelectLayerByLocation_management(CC_2010_2015, 'INTERSECT', lidar_2015,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2015, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2010_to_2015_lidar_2015.shp')

sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

## 2016
CC_5_years_before_2016 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2011 and "HARVYR" < 2016')
arcpy.CopyFeatures_management(CC_5_years_before_2016, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2011_2016.shp')

sweden_lidar_2016 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2016')
arcpy.CopyFeatures_management(sweden_lidar_2016, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2016.shp')

CC_2011_2016 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_2011_2016.shp'
lidar_2016 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2016.shp'

harvested_5_years_before_lidar_2016 = arcpy.SelectLayerByLocation_management(CC_2011_2016, 'INTERSECT', lidar_2016,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2016, 'Y:/William/RIP_NATURE_2/sweden/CC_Year/CC_2011_to_2016_lidar_2016.shp')

# Merge
inputs = 'Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2004_to_2009_lidar_2009.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2005_to_2010_lidar_2010.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2006_to_2011_lidar_2011.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2007_to_2012_lidar_2012shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2008_to_2013_lidar_2013.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2009_to_2014_lidar_2014.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2010_to_2015_lidar_2015.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2011_to_2016_lidar_2016.shp'
output = 'Y:/William/RIP_NATURE_2/sweden/sweden_harvest_before_LiDAR_merged.shp'
arcpy.management.Merge(inputs, output)
# Merge
inputs = 'Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2004_to_2009_lidar_2009.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2005_to_2010_lidar_2010.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2006_to_2011_lidar_2011.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2007_to_2012_lidar_2012.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2008_to_2013_lidar_2013.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2009_to_2014_lidar_2014.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2010_to_2015_lidar_2015.shp','Y:/William/RIP_NATURE_2/Sweden/CC_Year/CC_2011_to_2016_lidar_2016.shp'
output = 'Y:/William/RIP_NATURE_2/sweden/sweden_harvest_before_LiDAR_merged.shp'
arcpy.management.Merge(inputs, output)
sweden_post_harvest_year = pd.DataFrame.spatial.from_featureclass('Y:/William/RIP_NATURE_2/sweden/sweden_harvest_before_LiDAR_merged.shp')
sweden_post_harvest_year.head()<jupyter_output><empty_output><jupyter_text>We were only interested in harvest sites that were near a stream. We buffered the selected harvested sites with 30 m. These buffered harveste sites where then used to intersect with the topographical stream network<jupyter_code>arcpy.Buffer_analysis('Y:/William/RIP_NATURE_2/sweden/sweden_harvest_before_LiDAR_merged.shp', 'Y:/William/RIP_NATURE_2/Sweden/Buffer_CC_after_LiDAR/Buffered_post_harvest_scanne_CCs.shp', '30 Meters')
infeatures = (topographical_streams_sweden, 'Y:/William/RIP_NATURE_2/Sweden/Buffer_CC_after_LiDAR/Buffered_post_harvest_scanne_CCs.shp')
intersected_streams = 'Y:/William/RIP_NATURE_2/Sweden/Intersected_Streams.shp'
arcpy.Intersect_analysis(infeatures, intersected_streams)<jupyter_output><empty_output><jupyter_text>Select and export streams longer than 20m. Stream sections had to be at least  >20m because boreal canopy gaps can be naturally 0.1 ha, so we wanted to avoid these potentially natural gaps<jupyter_code>intersected_streams = 'Y:/William/RIP_NATURE_2/Sweden/Intersected_Streams.shp'

arcpy.AddGeometryAttributes_management(intersected_streams, 'LENGTH', 'METERS')
attribute_selected = arcpy.SelectLayerByAttribute_management(intersected_streams, 'NEW_SELECTION', 
                                        '"LENGTH" > 20')
arcpy.CopyFeatures_management(attribute_selected, 'Y:/William/RIP_NATURE_2/Sweden/stream_segments_longer_than_20m.shp')<jupyter_output><empty_output><jupyter_text>Next the stream segments were buffered on both sides by 20 m and dissolved seperatly.<jupyter_code># Streams
stream_segments_over_20m = 'Y:/William/RIP_NATURE_2/Sweden/stream_segments_longer_than_20m.shp'

# buffer outputs
buffer_left_side = 'Y:/William/RIP_NATURE_2/Sweden/Buffer_CC_after_LiDAR/Buffer_left_harvested_sites.shp'
buffer_right_side = 'Y:/William/RIP_NATURE_2/Sweden/Buffer_CC_after_LiDAR/Buffer_right_harvested_sites.shp'

buffer_distance = '20 Meters'
sideType_left = 'LEFT'
sideType_right = 'RIGHT'
endType = 'FLAT'
dissolveType = 'NONE'
dissolveField = 'NONE'


# Buffer Left
arcpy.Buffer_analysis(stream_segments_over_20m, buffer_left_side, buffer_distance, sideType_left, endType, 
                      dissolveType)

# Buffer Right
arcpy.Buffer_analysis(stream_segments_over_20m, buffer_right_side, buffer_distance, sideType_right, endType, 
                      dissolveType)
buffer_left_side = 'Y:/William/RIP_NATURE_2/Sweden/Buffer_CC_after_LiDAR/Buffer_left_harvested_sites.shp'
buffer_right_side = 'Y:/William/RIP_NATURE_2/Sweden/Buffer_CC_after_LiDAR/Buffer_right_harvested_sites.shp'

dissolve_left = 'Y:/William/RIP_NATURE_2/Sweden/dissolve/dissolve_left.shp'
dissolve_right = 'Y:/William/RIP_NATURE_2/Sweden/dissolve/dissolve_right.shp'

# Dissolve and copy left buffers
arcpy.Dissolve_management(buffer_left_side, dissolve_left,"", "", "SINGLE_PART", "UNSPLIT_LINES")

# Dissolve and copy Right buffers
arcpy.Dissolve_management(buffer_right_side, dissolve_right,"", "", "SINGLE_PART", "UNSPLIT_LINES")<jupyter_output><empty_output><jupyter_text>Select and merge buffer polygons within 10 m of a harvesting site that were harvested before being lidar scanned.<jupyter_code>dissolve_left = 'Y:/William/RIP_NATURE_2/Sweden/dissolve/dissolve_left.shp'
dissolve_right = 'Y:/William/RIP_NATURE_2/Sweden/dissolve/dissolve_right.shp'

harvestsites =  'Y:/William/RIP_NATURE_2/Sweden/CC_Before_LiDAR_within30m_of_stream/CC_BeforeLiDAR_wthin30mof_Stream.shp'

copy_left = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/within_10m_of_CC_left_buffer.shp'
copy_right = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/within_10m_of_CC_right_buffer.shp'

selection_type = 'NEW_SELECTION'

# Left
select_Left = arcpy.SelectLayerByLocation_management(dissolve_left, 'WITHIN_A_DISTANCE', harvestsites, '10 Meters', selection_type)
arcpy.CopyFeatures_management(select_Left, copy_left)
# Right
select_Right = arcpy.SelectLayerByLocation_management(dissolve_right, 'WITHIN_A_DISTANCE', harvestsites, '10 Meters', selection_type)
arcpy.CopyFeatures_management(select_Right, copy_right)

# Dissolved and selected buffers
Dissolved_left_side_within_10m_of_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/CC_within_10mof_left_buffer.shp'
Dissolved_right_right_within_10m_of_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/CC_within_10mof_right_buffer.shp'

# Merged buffers
merged_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/merged_buffer_within_10m_of_CC.shp'
arcpy.Merge_management([Dissolved_left_side_within_10m_of_harvest_sites, Dissolved_right_right_within_10m_of_harvest_sites], merged_buffer_within_10m_of_CC)<jupyter_output><empty_output><jupyter_text>Dissolve merged buffers<jupyter_code># Merged polygons
merged_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/merged_buffer_within_10m_of_CC.shp'
# Dissolved polygons
Dissolved_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/merged_buffer_within_10m_of_CC_dissolved.shp'

arcpy.Dissolve_management(merged_buffer_within_10m_of_CC, Dissolved_buffer_within_10m_of_CC,"", "", "SINGLE_PART", "DISSOLVE_LINES")<jupyter_output><empty_output><jupyter_text>Finally the Swedish forest mask was used to clip the buffer zones to make sure that only forest land were included<jupyter_code>Dissolved_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Sweden/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/merged_buffer_within_10m_of_CC_dissolved.shp'
Swedis_forest_land = 'Y:/William/RIP_NATURE_2/Sweden/ForestryLand.shp'
Post_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/Sweden/buffers_post_harvest.shp'
arcpy.Clip_analysis(Dissolved_buffer_within_10m_of_CC , Swedis_forest_land, Post_harvest_buffer_zones)
nr_sites_post_harvest = arcpy.GetCount_management(Post_harvest_buffer_zones)
print(nr_sites_post_harvest)<jupyter_output>70550<jupyter_text>Extract canopy cover statistics within the buffer zones using zonal statistcs as table<jupyter_code>import arcpy
from arcpy.sa import *
Post_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/Sweden/buffers_post_harvest.shp'

# Set the local variables
inZoneData = Post_harvest_buffer_zones
zoneField = 'FID'
inValueRaster = 'Y:/William/RIP_NATURE_2/Sweden/Tradhojdsraster/Mosaic_Hela_1_0.tif'
outTable = 'Y:/William/RIP_NATURE_2/Sweden/ZonalStatistics/sweden_post_harvest_buffers.dbf'

# Execute ZonalStatisticsAsTable
ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster, 
                                 outTable, "NODATA", "SUM", "")<jupyter_output><empty_output><jupyter_text>Pre-Harvest stream buffers  To relate the post harvest canopy densities I also need pre harvest canopy densities. To this end we selected sites that were harvested after being scanned with LiDAR. This means that they were in their unharvested state during the LiDAR scan.<jupyter_code>sweden_tiles_index = 'Y:/William/RIP_NATURE_2/Sweden/LiDAR/swedish_lidar_selected.shp'
sweden_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/CC/Swedish_CC.shp'

#The earliest scan was done 2009. 
#### 2009 ####
# CC after 2009
CC_after_2009 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2009')
arcpy.CopyFeatures_management(CC_after_2009, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2009.shp')
# LiDAR 2015
sweden_lidar_2009 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2009')
arcpy.CopyFeatures_management(sweden_lidar_2009, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2009')

CC_after_2009 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2009.shp'
sweden_lidar_2009 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2009.shp'

CC_after_lidar_2009 = arcpy.SelectLayerByLocation_management(CC_after_2009, 'INTERSECT', sweden_lidar_2009,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2009, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2009.shp')

# CC after 2010
CC_after_2010 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2010')
arcpy.CopyFeatures_management(CC_after_2010, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2010.shp')
# LiDAR 2010
sweden_lidar_2010 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2010')
arcpy.CopyFeatures_management(sweden_lidar_2010, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2010')

CC_after_2010 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2010.shp'
sweden_lidar_2010 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2010.shp'

CC_after_lidar_2010 = arcpy.SelectLayerByLocation_management(CC_after_2010, 'INTERSECT', sweden_lidar_2010,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2010, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2010.shp')

# CC after 2011
CC_after_2011 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2011')
arcpy.CopyFeatures_management(CC_after_2011, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2011.shp')
# LiDAR 2011
sweden_lidar_2011 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2011')
arcpy.CopyFeatures_management(sweden_lidar_2011, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2011')

CC_after_2011 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2011.shp'
sweden_lidar_2011 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2011.shp'

CC_after_lidar_2011 = arcpy.SelectLayerByLocation_management(CC_after_2011, 'INTERSECT', sweden_lidar_2011,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2011, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2011.shp')

# CC after 2012
CC_after_2012 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2012')
arcpy.CopyFeatures_management(CC_after_2012, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2012.shp')
# LiDAR 2012
sweden_lidar_2012 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2012')
arcpy.CopyFeatures_management(sweden_lidar_2012, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2012')

CC_after_2012 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2012.shp'
sweden_lidar_2012 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2012.shp'

CC_after_lidar_2012 = arcpy.SelectLayerByLocation_management(CC_after_2012, 'INTERSECT', sweden_lidar_2012,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2012, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2012.shp')

# CC after 2013
CC_after_2013 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2013')
arcpy.CopyFeatures_management(CC_after_2013, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2013.shp')
# LiDAR 2013
sweden_lidar_2013 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2013')
arcpy.CopyFeatures_management(sweden_lidar_2013, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2013')

CC_after_2013 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2013.shp'
sweden_lidar_2013 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2013.shp'

CC_after_lidar_2013 = arcpy.SelectLayerByLocation_management(CC_after_2013, 'INTERSECT', sweden_lidar_2013,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2013, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2013.shp')
# CC after 2014
CC_after_2014 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2014')
arcpy.CopyFeatures_management(CC_after_2014, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2014.shp')
# LiDAR 2014
sweden_lidar_2014 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2014')
arcpy.CopyFeatures_management(sweden_lidar_2014, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2014')

CC_after_2014 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2014.shp'
sweden_lidar_2014 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2014.shp'

CC_after_lidar_2014 = arcpy.SelectLayerByLocation_management(CC_after_2014, 'INTERSECT', sweden_lidar_2014,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2014, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2014.shp')

# CC after 2015
CC_after_2015 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2015')
arcpy.CopyFeatures_management(CC_after_2015, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2015.shp')
# LiDAR 2015
sweden_lidar_2015 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2015')
arcpy.CopyFeatures_management(sweden_lidar_2015, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2015')

CC_after_2015 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2015.shp'
sweden_lidar_2015 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2015.shp'

CC_after_lidar_2015 = arcpy.SelectLayerByLocation_management(CC_after_2015, 'INTERSECT', sweden_lidar_2015,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2015, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2015.shp')

# CC after 2016
CC_after_2016 = arcpy.SelectLayerByAttribute_management(sweden_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2016')
arcpy.CopyFeatures_management(CC_after_2016, 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2016.shp')
# LiDAR 2016
sweden_lidar_2016 = arcpy.SelectLayerByAttribute_management(sweden_tiles_index, 'NEW_SELECTION', '"final_year" = 2016')
arcpy.CopyFeatures_management(sweden_lidar_2016, 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2016')

CC_after_2016 = 'Y:/William/RIP_NATURE_2/sweden/Post_harvest/CC_after_2016.shp'
sweden_lidar_2016 = 'Y:/William/RIP_NATURE_2/sweden/LiDAR_Years/sweden_lidar_2016.shp'

CC_after_lidar_2016 = arcpy.SelectLayerByLocation_management(CC_after_2016, 'INTERSECT', sweden_lidar_2016,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2016, 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2016.shp')
# Merge reference sites
CC_after_lidar_2009 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2009.shp'
CC_after_lidar_2010 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2010.shp'
CC_after_lidar_2011 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2011.shp'
CC_after_lidar_2012 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2012.shp'
CC_after_lidar_2013 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2013.shp'
CC_after_lidar_2014 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2014.shp'
CC_after_lidar_2015 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2015.shp'
CC_after_lidar_2016 = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/CC_after_lidar_2016.shp'

merge_sweden_reference_sites = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/merge_sweden_reference_sites.shp'
arcpy.Merge_management([CC_after_lidar_2009,CC_after_lidar_2010,CC_after_lidar_2011,CC_after_lidar_2012,CC_after_lidar_2013,CC_after_lidar_2014,CC_after_lidar_2015, CC_after_lidar_2016], merge_sweden_reference_sites)<jupyter_output><empty_output><jupyter_text>Streams from the Swedish property map were intersected with  unharvested sites<jupyter_code>topographical_streams_sweden = 'Y:/William/RIP_NATURE_2/Sweden/hl_riks.shp'
Relevant_pre_harvest_sites = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/merge_sweden_reference_sites.shp'
Intersected_reference_streams = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/streams_intersect_CC.shp'

infeatures = (topographical_streams_sweden, Relevant_pre_harvest_sites)
arcpy.Intersect_analysis(infeatures, Intersected_reference_streams)
import pandas as pd
import matplotlib as mpl
from arcgis.features import SpatialDataFrame
sweden_pre_harvest = 'Y:/William/RIP_NATURE_2/sweden/Pre_harvest/merge_sweden_reference_sites.shp'
sweden_pre_harvest_year = pd.DataFrame.spatial.from_featureclass(sweden_pre_harvest)

import matplotlib.pyplot as plt
sweden_pre_harvest_year = sweden_pre_harvest_year['HARVYR'].value_counts().sort_index(ascending=True)


COLOR = 'grey'
mpl.rcParams['text.color'] = COLOR
mpl.rcParams['axes.labelcolor'] = COLOR
mpl.rcParams['xtick.color'] = COLOR
mpl.rcParams['ytick.color'] = COLOR

plt.rcParams.update({'font.size': 22})

fig = sweden_pre_harvest_year.plot(kind='bar', 
                        alpha=0.75,
                        figsize=(20,10),
                        color=COLOR,
                        rot='0')
plt.ylabel("Number of harvest sites")
plt.xlabel("Harvest year")

plt.show()<jupyter_output><empty_output><jupyter_text>Stream segments over 20 m were buffered by 20 m and intersected with unharvested sites<jupyter_code>Intersected_reference_streams = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/streams_intersect_CC.shp'
Intersected_reference_streams_over_20m = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/streams_intersect_CC_over20m.shp'

# Add lenght
arcpy.AddGeometryAttributes_management(Intersected_reference_streams, 'LENGTH', 'METERS')

# Select and copy segments over 20 m
attribute_selected_ref = arcpy.SelectLayerByAttribute_management(Intersected_reference_streams, 'NEW_SELECTION', '"LENGTH" > 20')
arcpy.CopyFeatures_management(attribute_selected_ref, Intersected_reference_streams_over_20m)

Intersected_reference_streams_over_20m = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/streams_intersect_CC_over20m.shp'
Intersected_reference_streams_over_20m_buffered_by_20m = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/streams_intersect_CC_over20m_buffered_20m.shp'

Relevant_pre_harvest_sites = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/CC_2011_to_2016_LiDAR_2009_or_2010.shp'
Intersect_buffered_stream_segments_with_unharvested_sites = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/Buffered_reference_streams_intersect_CC.shp'

# Buffer stream segments by 20 m
arcpy.Buffer_analysis(Intersected_reference_streams_over_20m, Intersected_reference_streams_over_20m_buffered_by_20m, '20 Meters')
# Intersect Buffered stream segments with unharvest sites
infeatures = (Intersected_reference_streams_over_20m_buffered_by_20m, Relevant_pre_harvest_sites)
arcpy.Intersect_analysis(infeatures, Intersect_buffered_stream_segments_with_unharvested_sites)<jupyter_output><empty_output><jupyter_text>Finally the reference buffers were dissolved based on harvest site ID.<jupyter_code>Intersect_buffered_stream_segments_with_unharvested_sites = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/Buffered_reference_streams_intersect_CC.shp'
Dissolve_reference_buffer_zones = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'
dissolveFields = ["FID_CC_201"]
arcpy.Dissolve_management(Intersect_buffered_stream_segments_with_unharvested_sites, Dissolve_reference_buffer_zones, dissolveFields, "", "SINGLE_PART", "DISSOLVE_LINES")<jupyter_output><empty_output><jupyter_text>Extract canopy cover within reference buffer polygons<jupyter_code>import arcpy
from arcpy.sa import *
Post_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'

# Set the local variables
inZoneData = Post_harvest_buffer_zones
zoneField = 'FID'
inValueRaster = 'Y:/William/RIP_NATURE_2/Sweden/Tradhojdsraster/Mosaic_Hela_1_0.tif'
outTable = 'Y:/William/RIP_NATURE_2/Sweden/ZonalStatistics/sweden_pre_harvest_buffers.dbf'

# Execute ZonalStatisticsAsTable
ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster, 
                                 outTable, "NODATA", "SUM", "")<jupyter_output><empty_output><jupyter_text>Sweden harvest sites and reference sites<jupyter_code># Harvestes sites 
Harvested = 'Y:/William/RIP_NATURE_2/Sweden/buffers_post_harvest.shp'
print(arcpy.GetCount_management(Harvested),'Harvested buffers')


# Reference sites
Reference = 'Y:/William/RIP_NATURE_2/Sweden/Reference_sites/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'
print(arcpy.GetCount_management(Reference),'Reference unharvested buffers')<jupyter_output>70550 Harvested buffers
21249 Reference unharvested buffers<jupyter_text>Finland  Post Harvest Post harvest means that the harvest sites were scanned with LiDAR within 5 years after being harvested. For example a site harvested between 2003 and 2008 would be in its harvested condition during a lidar scan that was conducted 2008. The shapefile 'Y:/William/RIP_NATURE_2/Finland/LIDAR/LiDARScanYear.shp' contains the column "final_year" that indicates when a sites was scanned with LiDAR. The file 'Y:/RIP_NATURE/Finland/Finland/NewclearcutsFromAntti/CC_declaration25112020.gdb/CC_declaration25112020.gdb/MKI_CC' contains harvest polygons and a column "year" that indicate which year the sites were harvested. Harvests between 2003 and 2018 were used.<jupyter_code>import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.parallelProcessingFactor = "100%"
finland_tiles_index = 'Y:/William/RIP_NATURE_2/Finland/LIDAR/LiDARScanYear.shp'
finland_harvest_sites = 'Y:/RIP_NATURE/Finland/Finland/NewclearcutsFromAntti/CC_declaration25112020.gdb/CC_declaration25112020.gdb/MKI_CC'

## 2008
CC_5_years_before_2008 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2003 and "year" < 2008')
arcpy.CopyFeatures_management(CC_5_years_before_2008, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2003_2008.shp')

finland_lidar_2008 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2008')
arcpy.CopyFeatures_management(finland_lidar_2008, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2008.shp')

CC_2003_2008 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2003_2008.shp'
lidar_2008 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2008.shp'

harvested_5_years_before_lidar_2008 = arcpy.SelectLayerByLocation_management(CC_2003_2008, 'INTERSECT', lidar_2008,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2008, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2003_to_2008_lidar_2008.shp')

## 2009
CC_5_years_before_2009 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2004 and "year" < 2009')
arcpy.CopyFeatures_management(CC_5_years_before_2009, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2004_2009.shp')

finland_lidar_2009 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2009')
arcpy.CopyFeatures_management(finland_lidar_2009, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2009.shp')

CC_2004_2009 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2004_2009.shp'
lidar_2009 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2009.shp'

harvested_5_years_before_lidar_2009 = arcpy.SelectLayerByLocation_management(CC_2004_2009, 'INTERSECT', lidar_2009,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2009, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2004_to_2009_lidar_2009.shp')

## 2010
CC_5_years_before_2010 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2005 and "year" < 2010')
arcpy.CopyFeatures_management(CC_5_years_before_2010, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2005_2010.shp')

finland_lidar_2010 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2010')
arcpy.CopyFeatures_management(finland_lidar_2010, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2010.shp')

CC_2005_2010 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2005_2010.shp'
lidar_2010 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2010.shp'

harvested_5_years_before_lidar_2010 = arcpy.SelectLayerByLocation_management(CC_2005_2010, 'INTERSECT', lidar_2010,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2010, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2005_to_2010_lidar_2010.shp')

## 2011
CC_5_years_before_2011 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2006 and "year" < 2011')
arcpy.CopyFeatures_management(CC_5_years_before_2011, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2006_2011.shp')

finland_lidar_2011 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2011')
arcpy.CopyFeatures_management(finland_lidar_2011, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2011.shp')

CC_2006_2011 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2006_2011.shp'
lidar_2011 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2011.shp'

harvested_5_years_before_lidar_2011 = arcpy.SelectLayerByLocation_management(CC_2006_2011, 'INTERSECT', lidar_2011,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2011, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2006_to_2011_lidar_2011.shp')

## 2012
CC_5_years_before_2012 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2007 and "year" < 2012')
arcpy.CopyFeatures_management(CC_5_years_before_2012, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2007_2012.shp')

finland_lidar_2012 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2012')
arcpy.CopyFeatures_management(finland_lidar_2012, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2012.shp')

CC_2007_2012 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2007_2012.shp'
lidar_2012 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2012.shp'

harvested_5_years_before_lidar_2012 = arcpy.SelectLayerByLocation_management(CC_2007_2012, 'INTERSECT', lidar_2012,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2012, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2007_to_2012_lidar_2012.shp')

## 2013
CC_5_years_before_2013 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2008 and "year" < 2013')
arcpy.CopyFeatures_management(CC_5_years_before_2013, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2008_2013.shp')

finland_lidar_2013 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2013')
arcpy.CopyFeatures_management(finland_lidar_2013, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2013.shp')

CC_2008_2013 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2008_2013.shp'
lidar_2013 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2013.shp'

harvested_5_years_before_lidar_2013 = arcpy.SelectLayerByLocation_management(CC_2008_2013, 'INTERSECT', lidar_2013,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2013, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2008_to_2013_lidar_2013.shp')

## 2014
CC_5_years_before_2014 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2009 and "year" < 2014')
arcpy.CopyFeatures_management(CC_5_years_before_2014, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2009_2014.shp')

finland_lidar_2014 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2014')
arcpy.CopyFeatures_management(finland_lidar_2014, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2014.shp')

CC_2009_2014 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2009_2014.shp'
lidar_2014 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2014.shp'

harvested_5_years_before_lidar_2014 = arcpy.SelectLayerByLocation_management(CC_2009_2014, 'INTERSECT', lidar_2014,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2014, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2009_to_2014_lidar_2014.shp')

## 2015
CC_5_years_before_2015 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2010 and "year" < 2015')
arcpy.CopyFeatures_management(CC_5_years_before_2015, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2010_2015.shp')

finland_lidar_2015 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2015')
arcpy.CopyFeatures_management(finland_lidar_2015, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2015.shp')

CC_2010_2015 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2010_2015.shp'
lidar_2015 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2015.shp'

harvested_5_years_before_lidar_2015 = arcpy.SelectLayerByLocation_management(CC_2010_2015, 'INTERSECT', lidar_2015,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2015, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2010_to_2015_lidar_2015.shp')

## 2016
CC_5_years_before_2016 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2011 and "year" < 2016')
arcpy.CopyFeatures_management(CC_5_years_before_2016, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2011_2016.shp')

finland_lidar_2016 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2016')
arcpy.CopyFeatures_management(finland_lidar_2016, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2016.shp')

CC_2011_2016 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2011_2016.shp'
lidar_2016 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2016.shp'

harvested_5_years_before_lidar_2016 = arcpy.SelectLayerByLocation_management(CC_2011_2016, 'INTERSECT', lidar_2016,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2016, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2011_to_2016_lidar_2016.shp')

## 2017
CC_5_years_before_2017 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2012 and "year" < 2017')
arcpy.CopyFeatures_management(CC_5_years_before_2017, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2012_2017.shp')

finland_lidar_2017 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2017')
arcpy.CopyFeatures_management(finland_lidar_2017, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2017.shp')

CC_2012_2017 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2012_2017.shp'
lidar_2017 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2017.shp'

harvested_5_years_before_lidar_2017 = arcpy.SelectLayerByLocation_management(CC_2012_2017, 'INTERSECT', lidar_2017,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2017, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2012_to_2017_lidar_2017.shp')

## 2018
CC_5_years_before_2018 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2013 and "year" < 2018')
arcpy.CopyFeatures_management(CC_5_years_before_2018, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2013_2018.shp')

finland_lidar_2018 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2018')
arcpy.CopyFeatures_management(finland_lidar_2018, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2018.shp')

CC_2013_2018 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_2013_2018.shp'
lidar_2018 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2018.shp'

harvested_5_years_before_lidar_2018 = arcpy.SelectLayerByLocation_management(CC_2013_2018, 'INTERSECT', lidar_2018,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2018, 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2013_to_2018_lidar_2018.shp')<jupyter_output><empty_output><jupyter_text>All selected harvest sites were merged into one shapefile for further processing.<jupyter_code># Merge
CC_2003_to_2008_lidar_2008 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2003_to_2008_lidar_2008.shp'
CC_2004_to_2009_lidar_2009 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2004_to_2009_lidar_2009.shp'
CC_2005_to_2010_lidar_2010 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2005_to_2010_lidar_2010.shp'
CC_2006_to_2011_lidar_2011 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2006_to_2011_lidar_2011.shp'
CC_2007_to_2012_lidar_2012 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2007_to_2012_lidar_2012.shp'
CC_2008_to_2013_lidar_2013 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2008_to_2013_lidar_2013.shp'
CC_2009_to_2014_lidar_2014 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2009_to_2014_lidar_2014.shp'
CC_2010_to_2015_lidar_2015 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2010_to_2015_lidar_2015.shp'
CC_2011_to_2016_lidar_2016 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2011_to_2016_lidar_2016.shp'
CC_2012_to_2017_lidar_2017 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2012_to_2017_lidar_2017.shp'
CC_2013_to_2018_lidar_2018 = 'Y:/William/RIP_NATURE_2/Finland/CC_Year/CC_2013_to_2018_lidar_2018.shp'

inputs = CC_2003_to_2008_lidar_2008, CC_2004_to_2009_lidar_2009, CC_2005_to_2010_lidar_2010,CC_2006_to_2011_lidar_2011,CC_2007_to_2012_lidar_2012,CC_2008_to_2013_lidar_2013,CC_2009_to_2014_lidar_2014,CC_2010_to_2015_lidar_2015,CC_2011_to_2016_lidar_2016,CC_2012_to_2017_lidar_2017,CC_2013_to_2018_lidar_2018

output = 'Y:/William/RIP_NATURE_2/finland/finland_harvest_before_LiDAR_merged.shp'
arcpy.management.Merge(inputs, output)<jupyter_output><empty_output><jupyter_text>The selected harvest sites were intersected with streams selected from this stream dataset: https://ckan.ymparisto.fi/dataset/%7BE15CD0B3-3934-4ABC-BF23-A8C622FA6A57%7DThe stream were selected based on the upper catchment area of stream endpoints. "ValuLuokka" 1 and 2 represents small streams.<jupyter_code>Original_Streams_Finland = 'Y:/William/RIP_NATURE_2/Finland/streams/Uoma10.shp'
selected_streams_Finland =  'Y:/William/RIP_NATURE_2/Finland/streams/Uoma1and2.shp'

attribute_selected_ref = arcpy.SelectLayerByAttribute_management(Original_Streams_Finland, 'NEW_SELECTION', '"ValuLuokka" <= 2')
arcpy.CopyFeatures_management(attribute_selected_ref, selected_streams_Finland)<jupyter_output><empty_output><jupyter_text>harvest polygons within 30 m of a stream were selected and buffered by 30 m. These buffered harvest sites will be used later to intersect stream segments.<jupyter_code>import arcpy
Post_harvest_sites_Finland = 'Y:/William/RIP_NATURE_2/finland/finland_harvest_before_LiDAR_merged.shp'
selected_streams_Finland =  'Y:/William/RIP_NATURE_2/Finland/streams/Uoma1and2.shp'
CC_BeforeLiDAR_wthin30mof_Stream = 'Y:/William/RIP_NATURE_2/Finland/CCWithin5YearsBBefore_LiDAR_within30m_of_stream.shp'
Buffered_CC_BeforeLiDAR_wthin30mof_Stream = 'Y:/William/RIP_NATURE_2/Finland/BufferedCCWithin5YearsBBefore_LiDAR_within30m_of_stream.shp'

# select harvest sites within 30 m of streams
selected_within_30m_stream_Finland = arcpy.SelectLayerByLocation_management(Post_harvest_sites_Finland, 'WITHIN_A_DISTANCE', selected_streams_Finland, '30 Meters')
# Export selected harvest sites to new shapefile
arcpy.CopyFeatures_management(selected_within_30m_stream_Finland, CC_BeforeLiDAR_wthin30mof_Stream)

# Buffer harvest sites within 30 m of streams by 30 m
arcpy.Buffer_analysis(CC_BeforeLiDAR_wthin30mof_Stream, Buffered_CC_BeforeLiDAR_wthin30mof_Stream, '30 Meters')<jupyter_output><empty_output><jupyter_text>The buffered harvest sites were intersected with the selected streams. Intersected stream segments that were longer than 20 m were selected for further processing.<jupyter_code># Inputs for intersect
selected_streams_Finland =  'Y:/William/RIP_NATURE_2/Finland/streams/Uoma1and2.shp'
Buffered_CC_BeforeLiDAR_wthin30mof_Stream = 'Y:/William/RIP_NATURE_2/Finland/BufferedCCWithin5YearsBBefore_LiDAR_within30m_of_stream.shp'
intersected_streams_Finland = 'Y:/William/RIP_NATURE_2/Finland/Intersected_Streams.shp'

infeatures = (selected_streams_Finland, Buffered_CC_BeforeLiDAR_wthin30mof_Stream)
arcpy.Intersect_analysis(infeatures, intersected_streams_Finland)

# Add lenght to the stream segments
arcpy.AddGeometryAttributes_management(intersected_streams_Finland, 'LENGTH', 'METERS')

# Select segments longer than 20 m
attribute_selected = arcpy.SelectLayerByAttribute_management(intersected_streams_Finland, 'NEW_SELECTION', 
                                        '"LENGTH" > 20')
stream_segments_over_20m = 'Y:/William/RIP_NATURE_2/Finland/Intersected_Streams_longer_than_20m.shp'
arcpy.CopyFeatures_management(attribute_selected, stream_segments_over_20m)<jupyter_output><empty_output><jupyter_text>Next the stream segments where buffered and dissolved seperatly<jupyter_code># Streams
stream_segments_over_20m = 'Y:/William/RIP_NATURE_2/Finland/Intersected_Streams_longer_than_20m.shp'

# buffer outputs
buffer_left_side = 'Y:/William/RIP_NATURE_2/Finland/Buffer_CC_after_LiDAR/Buffer_left_harvested_sites.shp'
buffer_right_side = 'Y:/William/RIP_NATURE_2/Finland/Buffer_CC_after_LiDAR/Buffer_right_harvested_sites.shp'

buffer_distance = '20 Meters'
sideType_left = 'LEFT'
sideType_right = 'RIGHT'
endType = 'FLAT'
dissolveType = 'NONE'
dissolveField = 'NONE'

# Buffer Left and right sides seperatly
arcpy.Buffer_analysis(stream_segments_over_20m, buffer_left_side, buffer_distance, sideType_left, endType, dissolveType)
arcpy.Buffer_analysis(stream_segments_over_20m, buffer_right_side, buffer_distance, sideType_right, endType, dissolveType)

# DISOLVE
dissolve_left = 'Y:/William/RIP_NATURE_2/Finland/dissolve/dissolve_left.shp'
dissolve_right = 'Y:/William/RIP_NATURE_2/Finland/dissolve/dissolve_right.shp'

# Dissolve and copy both sides
arcpy.Dissolve_management(buffer_left_side, dissolve_left,"", "", "SINGLE_PART", "UNSPLIT_LINES")
arcpy.Dissolve_management(buffer_right_side, dissolve_right,"", "", "SINGLE_PART", "UNSPLIT_LINES")<jupyter_output><empty_output><jupyter_text>Select left and right buffers within 10 m from the unbuffered harvest sites that were located within 30m of a stream. Merge selected buffers into a single file.<jupyter_code>dissolve_left = 'Y:/William/RIP_NATURE_2/Finland/dissolve/dissolve_left.shp'
dissolve_right = 'Y:/William/RIP_NATURE_2/Finland/dissolve/dissolve_right.shp'

# Within 30 m of a stream
Post_Harvest_Finland = 'Y:/William/RIP_NATURE_2/Finland/CCWithin5YearsBBefore_LiDAR_within30m_of_stream.shp'

copy_left = 'Y:/William/RIP_NATURE_2/Finland/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/within_10m_of_CC_left_buffer.shp'
copy_right = 'Y:/William/RIP_NATURE_2/Finland/Dissolved_buffers_within_10m_of_stream_lidar_after_harvest/within_10m_of_CC_right_buffer.shp'

selection_type = 'NEW_SELECTION'

# Left
select_Left = arcpy.SelectLayerByLocation_management(dissolve_left, 'WITHIN_A_DISTANCE', Post_Harvest_Finland, '10 Meters', selection_type)
arcpy.CopyFeatures_management(select_Left, copy_left)
# Right
select_Right = arcpy.SelectLayerByLocation_management(dissolve_right, 'WITHIN_A_DISTANCE', Post_Harvest_Finland, '10 Meters', selection_type)
arcpy.CopyFeatures_management(select_Right, copy_right)


# Merged buffers
merged_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Finland/merged_buffer_within_10m_of_CC.shp'
arcpy.Merge_management([copy_left, copy_right], merged_buffer_within_10m_of_CC)<jupyter_output><empty_output><jupyter_text>Dissolve the merged buffers and clip with finish forest land to avoid buffers on aggricultural land.<jupyter_code># Merged polygons
merged_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Finland/merged_buffer_within_10m_of_CC.shp'
# Dissolved polygons
Dissolved_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Finland/merged_buffer_within_10m_of_CC_dissolved.shp'

arcpy.Dissolve_management(merged_buffer_within_10m_of_CC, Dissolved_buffer_within_10m_of_CC,"", "", "SINGLE_PART", "DISSOLVE_LINES")


finish_forest_land = 'Y:/William/RIP_NATURE_2/Finland/finish_forestry_mask_v2.gdb/finish_forestry_mask_v2'
Post_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/Finland/buffers_post_harvest.shp'
arcpy.Clip_analysis(Dissolved_buffer_within_10m_of_CC , finish_forest_land, Post_harvest_buffer_zones)<jupyter_output><empty_output><jupyter_text>Pre Harvest Pre harvest means that these sites were harvested after being scanned with LiDAR. They were thus intact during the liDAR scan and will be used as reference. Further we are only interested in buffers that intersect harvest sites. Additionally the buffers should be from stream segments longer than 20 m.<jupyter_code>arcpy.env.overwriteOutput = True
finland_tiles_index = 'Y:/William/RIP_NATURE_2/Finland/LIDAR/LiDARScanYear.shp'
finland_harvest_sites = 'Y:/William/RIP_NATURE_2/Finland/NewclearcutsFromAntti/CC_declaration25112020.gdb/CC_declaration25112020.gdb/MKI_CC'

#### 2008 ####
# CC 2008
CC_after_2008 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2008')
arcpy.CopyFeatures_management(CC_after_2008, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2008.shp')
# LiDAR 20108
finland_lidar_2008 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2008')
arcpy.CopyFeatures_management(finland_lidar_2008, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2008')

CC_after_2008 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2008.shp'
finland_lidar_2008 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2008.shp'

CC_after_lidar_2008 = arcpy.SelectLayerByLocation_management(CC_after_2008, 'INTERSECT', finland_lidar_2008,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2008, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2008.shp')

#### 2009 ####
# CC 2009
CC_after_2009 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2009')
arcpy.CopyFeatures_management(CC_after_2009, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2009.shp')
# LiDAR 2009
finland_lidar_2009 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2009')
arcpy.CopyFeatures_management(finland_lidar_2009, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2009')

CC_after_2009 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2009.shp'
finland_lidar_2009 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2009.shp'

CC_after_lidar_2009 = arcpy.SelectLayerByLocation_management(CC_after_2009, 'INTERSECT', finland_lidar_2009,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2009, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2009.shp')

#### 2010 ####
# CC 2010
CC_after_2010 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2010')
arcpy.CopyFeatures_management(CC_after_2010, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2010.shp')
# LiDAR 2010
finland_lidar_2010 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2010')
arcpy.CopyFeatures_management(finland_lidar_2010, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2010')

CC_after_2010 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2010.shp'
finland_lidar_2010 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2010.shp'

CC_after_lidar_2010 = arcpy.SelectLayerByLocation_management(CC_after_2010, 'INTERSECT', finland_lidar_2010,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2010, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2010.shp')

#### 2011 ####
# CC 2011
CC_after_2011 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2011')
arcpy.CopyFeatures_management(CC_after_2011, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2011.shp')
# LiDAR 2011
finland_lidar_2011 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2011')
arcpy.CopyFeatures_management(finland_lidar_2011, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2011')

CC_after_2011 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2011.shp'
finland_lidar_2011 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2011.shp'

CC_after_lidar_2011 = arcpy.SelectLayerByLocation_management(CC_after_2011, 'INTERSECT', finland_lidar_2011,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2011, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2011.shp')

#### 2012 ####
# CC 2012
CC_after_2012 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2012')
arcpy.CopyFeatures_management(CC_after_2012, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2012.shp')
# LiDAR 2012
finland_lidar_2012 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2012')
arcpy.CopyFeatures_management(finland_lidar_2012, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2012')

CC_after_2012 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2012.shp'
finland_lidar_2012 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2012.shp'

CC_after_lidar_2012 = arcpy.SelectLayerByLocation_management(CC_after_2012, 'INTERSECT', finland_lidar_2012,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2012, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2012.shp')

#### 2013 ####
# CC 2013
CC_after_2013 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2013')
arcpy.CopyFeatures_management(CC_after_2013, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2013.shp')
# LiDAR 2013
finland_lidar_2013 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2013')
arcpy.CopyFeatures_management(finland_lidar_2013, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2013')

CC_after_2013 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2013.shp'
finland_lidar_2013 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2013.shp'

CC_after_lidar_2013 = arcpy.SelectLayerByLocation_management(CC_after_2013, 'INTERSECT', finland_lidar_2013,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2013, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2013.shp')

#### 2014 ####
# CC 2014
CC_after_2014 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2014')
arcpy.CopyFeatures_management(CC_after_2014, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2014.shp')
# LiDAR 2014
finland_lidar_2014 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2014')
arcpy.CopyFeatures_management(finland_lidar_2014, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2014')

CC_after_2014 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2014.shp'
finland_lidar_2014 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2014.shp'

CC_after_lidar_2014 = arcpy.SelectLayerByLocation_management(CC_after_2014, 'INTERSECT', finland_lidar_2014,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2014, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2014.shp')

arcpy.env.overwriteOutput = True
finland_tiles_index = 'Y:/William/RIP_NATURE_2/Finland/LIDAR/LiDARScanYear.shp'
finland_harvest_sites = 'Y:/William/RIP_NATURE_2/Finland/NewclearcutsFromAntti/CC_declaration25112020.gdb/CC_declaration25112020.gdb/MKI_CC'


#### 2015 ####
# CC 2015
CC_after_2015 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2015')
arcpy.CopyFeatures_management(CC_after_2015, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2015.shp')
# LiDAR 2015
finland_lidar_2015 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2015')
arcpy.CopyFeatures_management(finland_lidar_2015, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2015')

CC_after_2015 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2015.shp'
finland_lidar_2015 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2015.shp'

CC_after_lidar_2015 = arcpy.SelectLayerByLocation_management(CC_after_2015, 'INTERSECT', finland_lidar_2015,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2015, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2015.shp')


#### 2016 ####
# CC 2016
CC_after_2016 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2016')
arcpy.CopyFeatures_management(CC_after_2016, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2016.shp')
# LiDAR 2016
finland_lidar_2016 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2016')
arcpy.CopyFeatures_management(finland_lidar_2016, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2016')

CC_after_2016 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2016.shp'
finland_lidar_2016 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2016.shp'

CC_after_lidar_2016 = arcpy.SelectLayerByLocation_management(CC_after_2016, 'INTERSECT', finland_lidar_2016,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2016, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2016.shp')


#### 2017 ####
# CC 2017
CC_after_2017 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2017')
arcpy.CopyFeatures_management(CC_after_2017, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2017.shp')
# LiDAR 2017
finland_lidar_2017 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2017')
arcpy.CopyFeatures_management(finland_lidar_2017, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2017')

CC_after_2017 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2017.shp'
finland_lidar_2017 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2017.shp'

CC_after_lidar_2017 = arcpy.SelectLayerByLocation_management(CC_after_2017, 'INTERSECT', finland_lidar_2017,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2017, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2017.shp')



#### 2018 ####
# CC 2018
CC_after_2018 = arcpy.SelectLayerByAttribute_management(finland_harvest_sites, 'NEW_SELECTION', '"year" > 2018')
arcpy.CopyFeatures_management(CC_after_2018, 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2018.shp')
# LiDAR 2018
finland_lidar_2018 = arcpy.SelectLayerByAttribute_management(finland_tiles_index, 'NEW_SELECTION', '"final_year" = 2018')
arcpy.CopyFeatures_management(finland_lidar_2018, 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2018')

CC_after_2018 = 'Y:/William/RIP_NATURE_2/finland/Post_harvest/CC_after_2018.shp'
finland_lidar_2018 = 'Y:/William/RIP_NATURE_2/finland/LiDAR_Years/finland_lidar_2018.shp'

CC_after_lidar_2018 = arcpy.SelectLayerByLocation_management(CC_after_2018, 'INTERSECT', finland_lidar_2018,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2018, 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2018.shp')



CC_after_lidar_2008 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2008.shp'
CC_after_lidar_2009 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2009.shp'
CC_after_lidar_2010 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2010.shp'
CC_after_lidar_2011 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2011.shp'
CC_after_lidar_2012 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2012.shp'
CC_after_lidar_2013 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2013.shp'
CC_after_lidar_2014 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2014.shp'
CC_after_lidar_2015 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2015.shp'
CC_after_lidar_2016 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2016.shp'
CC_after_lidar_2017 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2017.shp'
CC_after_lidar_2018 = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/CC_after_lidar_2018.shp'

inputs = CC_after_lidar_2008,CC_after_lidar_2009,CC_after_lidar_2010,CC_after_lidar_2011, CC_after_lidar_2012, CC_after_lidar_2013,CC_after_lidar_2014, CC_after_lidar_2015, CC_after_lidar_2016, CC_after_lidar_2017, CC_after_lidar_2018
output = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/finland_harvest_after_LiDAR_merged.shp'
arcpy.management.Merge(inputs, output)
# inputs
Pre_Harvest_Finland = 'Y:/William/RIP_NATURE_2/finland/pre_harvest/finland_harvest_after_LiDAR_merged.shp'
selected_streams_Finland =  'Y:/William/RIP_NATURE_2/Finland/streams/Uoma1and2.shp'

# output
Intersected_reference_streams_Finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Intersected_reference_streams.shp'

# Intersect
infeatures = (selected_streams_Finland , Pre_Harvest_Finland)
arcpy.Intersect_analysis(infeatures, Intersected_reference_streams_Finland)

# add geometry and select stream segments longer than 20 m.
Intersected_reference_streams_Finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Intersected_reference_streams.shp'
Intersected_reference_streams_over_20m_Finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Intersected_reference_streams_over_20m.shp'
# Add lenght
arcpy.AddGeometryAttributes_management(Intersected_reference_streams_Finland, 'LENGTH', 'METERS')

# Select and copy segments over 20 m
attribute_selected_ref = arcpy.SelectLayerByAttribute_management(Intersected_reference_streams_Finland, 'NEW_SELECTION', '"LENGTH" > 20')
arcpy.CopyFeatures_management(attribute_selected_ref, Intersected_reference_streams_over_20m_Finland)<jupyter_output><empty_output><jupyter_text>The 20 m segments were buffered by 20 m and intersected with and dissolved by the unharvested sites.<jupyter_code>import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.parallelProcessingFactor = "100%"


Intersected_reference_streams_over_20m_Finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Intersected_reference_streams_over_20m.shp'
Intersected_reference_streams_over_20m_Finland_Buffered = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Intersected_reference_streams_over_20m_buffered_by_20m.shp'


# Buffer stream segments by 20 m
arcpy.Buffer_analysis(Intersected_reference_streams_over_20m_Finland, Intersected_reference_streams_over_20m_Finland_Buffered, '20 Meters')

Pre_Harvest_Finland = 'Y:/William/RIP_NATURE_2/Finland/CC_After_LiDAR/CCAfterLiDAR.shp'
Intersect_buffered_stream_segments_with_unharvested_sites_Finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Buffered_reference_streams_intersect_CC.shp'


# Intersect Buffered stream segments with unharvest sites
infeatures = (Intersected_reference_streams_over_20m_Finland_Buffered, Pre_Harvest_Finland)
arcpy.Intersect_analysis(infeatures, Intersect_buffered_stream_segments_with_unharvested_sites_Finland)

# DISSOLVE BY HARVEST ID
Dissolve_reference_buffer_zones_Finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'
dissolveFields = ["OBJECTID_1"]
arcpy.Dissolve_management(Intersect_buffered_stream_segments_with_unharvested_sites_Finland, Dissolve_reference_buffer_zones_Finland, dissolveFields, "", "SINGLE_PART", "DISSOLVE_LINES")
# Harvestes sites 
Harvested_finland = 'Y:/William/RIP_NATURE_2/Finland/buffers_post_harvest.shp'
print(arcpy.GetCount_management(Harvested_finland),'Harvested buffers')


# Reference sites
Reference_finland = 'Y:/William/RIP_NATURE_2/Finland/Reference_sites/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'
print(arcpy.GetCount_management(Reference_finland),'Reference unharvested buffers')<jupyter_output>5795 Harvested buffers
2459 Reference unharvested buffers<jupyter_text>Copy finish canopy files to single folder<jupyter_code>import arcpy
import os
import glob
single_folder = 'Y:/NewCanopyFromFinland/Finland/CHM_single_folder/'
finish_canopy = glob.glob('Y:/NewCanopyFromFinland/Finland/CHM/**/*.tif', recursive = True)
print(len(finish_canopy))
for file in finish_canopy:
    head, tail = os.path.split(file)
    output = single_folder + tail
    arcpy.CopyRaster_management(file, output)<jupyter_output>27567<jupyter_text>Canada  Post Harvest<jupyter_code># number of reference harvest sites
nr_sites_post_harvest = arcpy.GetCount_management('Y:/William/RIP_NATURE_2/Canada/Harvestfrom2012.shp')
print('Number of harvests in NB between 2012 and 2018:',nr_sites_post_harvest)<jupyter_output>Number of harvests in NB between 2012 and 2018: 20271<jupyter_text>Select sites that were harvested within 5 years before being scanned with LiDARA shapeile with the lidar tiles was downloaded from http://www.snb.ca/geonb1/e/dc/catalogue-E.asp and is stored at: "Y:\William\RIP_NATURE_2\Canada\LiDAR_Years\geonb_li-idl_shp\geonb_li_idl_cgvd2013.shp". note that some areas where scanned both in 2017 and 2018. 2018 were used as final year and compared to harvests 2019 and onwards. for the reference sites 2016 and older were used. in the yeanr field they are set to 2018.<jupyter_code>nb_tiles_index = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/geonb_li-idl_shp/geonb_li_idl_cgvd2013.shp'
nb_harvest_sites = 'Y:/William/RIP_NATURE_2/Canada/Harvestfrom2012.shp'

# 2015
CC_5_years_before_2015 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2010 and "HARVYR" < 2015')
arcpy.CopyFeatures_management(CC_5_years_before_2015, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2010_2015.shp')

nb_lidar_2015 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2015')
arcpy.CopyFeatures_management(nb_lidar_2015, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2015')

CC_2010_2015 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2010_2015.shp'
lidar_2015 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2015.shp'

harvested_5_years_before_lidar_2015 = arcpy.SelectLayerByLocation_management(CC_2010_2015, 'INTERSECT', lidar_2015,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2015, 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2010_to_2015_lidar_2015.shp')

# 2016
CC_5_years_before_2016 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2011 and "HARVYR" < 2016')
arcpy.CopyFeatures_management(CC_5_years_before_2016, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2011_2016.shp')

nb_lidar_2016 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2016')
arcpy.CopyFeatures_management(nb_lidar_2016, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2016')

CC_2011_2016 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2011_2016.shp'
lidar_2016 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2016.shp'

harvested_5_years_before_lidar_2016 = arcpy.SelectLayerByLocation_management(CC_2011_2016, 'INTERSECT', lidar_2016,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2016, 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2011_to_2016_lidar_2016.shp')

# 2017
CC_5_years_before_2017 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2012 and "HARVYR" < 2017')
arcpy.CopyFeatures_management(CC_5_years_before_2017, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2012_2017.shp')

nb_lidar_2017 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2017')
arcpy.CopyFeatures_management(nb_lidar_2017, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2017')

CC_2012_2017 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2012_2017.shp'
lidar_2017 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2017.shp'

harvested_5_years_before_lidar_2017 = arcpy.SelectLayerByLocation_management(CC_2012_2017, 'INTERSECT', lidar_2017,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2017, 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2012_to_2017_lidar_2017.shp')

# 2018
CC_5_years_before_2018 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2012 and "HARVYR" < 2018')
arcpy.CopyFeatures_management(CC_5_years_before_2018, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2012_2018.shp')

nb_lidar_2018 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2018')
arcpy.CopyFeatures_management(nb_lidar_2018, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2018')

CC_2012_2018 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2012_2018.shp'
lidar_2018 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2018.shp'

harvested_5_years_before_lidar_2018 = arcpy.SelectLayerByLocation_management(CC_2012_2018, 'INTERSECT', lidar_2018,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(harvested_5_years_before_lidar_2018, 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2012_to_2018_lidar_2018.shp')

# Merge
inputs = 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2010_to_2015_lidar_2015.shp', 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2011_to_2016_lidar_2016.shp', 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2012_to_2017_lidar_2017.shp', 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_2012_to_2018_lidar_2018.shp'
output = 'Y:/William/RIP_NATURE_2/Canada/NB_harvest_before_LiDAR_merged.shp'
arcpy.management.Merge(inputs, output)<jupyter_output><empty_output><jupyter_text>Select sites within 30 m of a stream and buffer them with 30 m<jupyter_code>sites_scanned_after_harvest = 'Y:/William/RIP_NATURE_2/Canada/NB_harvest_before_LiDAR_merged.shp'
streams = 'Y:/William/RIP_NATURE_2/Canada/streams/PropertyStreams.shp'

#sites_near_streams = arcpy.SelectLayerByLocation_management(sites_scanned_after_harvest, 'WITHIN_A_DISTANCE', streams,'30 Meter', ) 
#arcpy.CopyFeatures_management(sites_near_streams, 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_merge_within_30m_of_stream.shp')

# Buffer harvest sites within 30 m of streams by 30 m
nb_cc_before_lidar_within_30m_of_stream = 'Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_merge_within_30m_of_stream.shp'
nb_buffered_CC_before_LiDAR_wthin_30m_of_Stream = 'Y:/William/RIP_NATURE_2/Canada/buffered_cc_within_5_years_before_LiDAR_within_30m_of_stream.shp'

arcpy.Buffer_analysis(nb_cc_before_lidar_within_30m_of_stream, nb_buffered_CC_before_LiDAR_wthin_30m_of_Stream, '30 Meters')

# number of reference harvest sites
nr_sites_lidar_post_harvest_near_stream = arcpy.GetCount_management('Y:/William/RIP_NATURE_2/Canada/CC_Year/CC_merge_within_30m_of_stream.shp')
print('Number of harvests in NB between that were scanned with LiDAR after harvest and located within 30m of a stream:', nr_sites_lidar_post_harvest_near_stream)<jupyter_output>Number of harvests in NB between that were scanned with LiDAR after harvest and located within 30m of a stream: 4044<jupyter_text>Stream segments intersected with buffered harvest sites that were longer than 20 m were selected for further processing<jupyter_code># Inputs for intersect
streams = 'Y:/William/RIP_NATURE_2/Canada/streams/PropertyStreams.shp'
nb_buffered_CC_before_LiDAR_wthin_30m_of_Stream = 'Y:/William/RIP_NATURE_2/Canada/buffered_cc_within_5_years_before_LiDAR_within_30m_of_stream.shp'
nb_intersected_streams = 'Y:/William/RIP_NATURE_2/Canada/Intersected_Streams.shp'

infeatures = (streams, nb_buffered_CC_before_LiDAR_wthin_30m_of_Stream)
arcpy.Intersect_analysis(infeatures, nb_intersected_streams)

# Add lenght to the stream segments
arcpy.AddGeometryAttributes_management(nb_intersected_streams, 'LENGTH', 'METERS')

# Select segments longer than 20 m
attribute_selected = arcpy.SelectLayerByAttribute_management(nb_intersected_streams, 'NEW_SELECTION', 
                                        '"LENGTH" > 20')
stream_segments_over_20m = 'Y:/William/RIP_NATURE_2/Canada/Intersected_Streams_longer_than_20m.shp'
arcpy.CopyFeatures_management(attribute_selected, stream_segments_over_20m)<jupyter_output><empty_output><jupyter_text>The stream segments were buffered and dissolved seperetly<jupyter_code># Streams
stream_segments_over_20m = 'Y:/William/RIP_NATURE_2/canada/Intersected_Streams_longer_than_20m.shp'


# buffer outputs
buffer_left_side = 'Y:/William/RIP_NATURE_2/Canada/Buffer_CC_after_LiDAR/Buffer_left_harvested_sites.shp'
buffer_right_side = 'Y:/William/RIP_NATURE_2/Canada/Buffer_CC_after_LiDAR/Buffer_right_harvested_sites.shp'

buffer_distance = '20 Meters'
sideType_left = 'LEFT'
sideType_right = 'RIGHT'
endType = 'FLAT'
dissolveType = 'NONE'
dissolveField = 'NONE'


# Buffer Left
arcpy.Buffer_analysis(stream_segments_over_20m, buffer_left_side, buffer_distance, sideType_left, endType, 
                      dissolveType)

# Buffer Right
arcpy.Buffer_analysis(stream_segments_over_20m, buffer_right_side, buffer_distance, sideType_right, endType, 
                      dissolveType)

# DISOLVE
dissolve_left = 'Y:/William/RIP_NATURE_2/Canada/dissolve/dissolve_left.shp'
dissolve_right = 'Y:/William/RIP_NATURE_2/Canada/dissolve/dissolve_right.shp'

# Dissolve and copy left buffers
arcpy.Dissolve_management(buffer_left_side, dissolve_left,"", "", "SINGLE_PART", "UNSPLIT_LINES")

# Dissolve and copy Right buffers
arcpy.Dissolve_management(buffer_right_side, dissolve_right,"", "", "SINGLE_PART", "UNSPLIT_LINES")<jupyter_output><empty_output><jupyter_text>Select left and right buffers within 10 m of harvest sites and merge them<jupyter_code>dissolve_left = 'Y:/William/RIP_NATURE_2/canada/dissolve/dissolve_left.shp'
dissolve_right = 'Y:/William/RIP_NATURE_2/canada/dissolve/dissolve_right.shp'


Post_Harvest_canada = 'Y:/William/RIP_NATURE_2/Canada/NB_harvest_before_LiDAR_merged.shp'

copy_left = 'Y:/William/RIP_NATURE_2/Canada/Dissolve/within_10m_of_CC_left_buffer.shp'
copy_right = 'Y:/William/RIP_NATURE_2/Canada/Dissolve/within_10m_of_CC_right_buffer.shp'

selection_type = 'NEW_SELECTION'

# Left
select_Left = arcpy.SelectLayerByLocation_management(dissolve_left, 'WITHIN_A_DISTANCE', Post_Harvest_canada, '10 Meters', selection_type)
arcpy.CopyFeatures_management(select_Left, copy_left)
# Right
select_Right = arcpy.SelectLayerByLocation_management(dissolve_right, 'WITHIN_A_DISTANCE', Post_Harvest_canada, '10 Meters', selection_type)
arcpy.CopyFeatures_management(select_Right, copy_right)


# Merged buffers
merged_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Canada/dissolve/merged_buffer_within_10m_of_CC.shp'
arcpy.Merge_management([copy_left, copy_right], merged_buffer_within_10m_of_CC)<jupyter_output><empty_output><jupyter_text>Dissolve the merged buffers and clip with harvest land<jupyter_code># Merged polygons
merged_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Canada/dissolve/merged_buffer_within_10m_of_CC.shp'

# Dissolved polygons
Dissolved_buffer_within_10m_of_CC = 'Y:/William/RIP_NATURE_2/Canada/dissolve/merged_buffer_within_10m_of_CC_dissolved.shp'

arcpy.Dissolve_management(merged_buffer_within_10m_of_CC, Dissolved_buffer_within_10m_of_CC,"", "", "SINGLE_PART", "DISSOLVE_LINES")

nb_forest_land = 'Y:/Canada/geonb_forest-foret_shp/geonb_forest-foret.shp'
Lidar_Post_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/canada/nb_buffers_lidar_post_harvest.shp'
arcpy.Clip_analysis(Dissolved_buffer_within_10m_of_CC , nb_forest_land, Lidar_Post_harvest_buffer_zones)<jupyter_output><empty_output><jupyter_text>Pre Harvest To evaluate canopy density after harvesting we need a estimate of unharvested canopy densities. To make the comparison as fair as possible we seleted sites that were harvested after being scanned with LiDAR. This means that the sites will be in unharvested conditions in the lidardata.Example: A site that was scanned with LiDAR 2016 and harvested 2017 would be intact in the lidar data.<jupyter_code>nb_tiles_index = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/geonb_li-idl_shp/geonb_li_idl_cgvd2013.shp'
nb_harvest_sites = 'Y:/William/RIP_NATURE_2/Canada/Harvestfrom2012.shp'

#### 2015 ####
# CC 2015
CC_after_2015 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2015')
arcpy.CopyFeatures_management(CC_after_2015, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_2015.shp')
# LiDAR 2015
nb_lidar_2015 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2015')
arcpy.CopyFeatures_management(nb_lidar_2015, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2015')

CC_after_2015 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_after_2015.shp'
nb_lidar_2015 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2015.shp'

CC_after_lidar_2015 = arcpy.SelectLayerByLocation_management(CC_after_2015, 'INTERSECT', nb_lidar_2015,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2015, 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/CC_after_lidar_2015.shp')

#### 2016 ####
# CC 2016
CC_after_2016 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2016')
arcpy.CopyFeatures_management(CC_after_2016, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_after_2016.shp')
# LiDAR 2016
nb_lidar_2016 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2016')
arcpy.CopyFeatures_management(nb_lidar_2016, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2016.shp')

CC_after_2016 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_after_2016.shp'
nb_lidar_2016 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2016.shp'

CC_after_lidar_2016 = arcpy.SelectLayerByLocation_management(CC_after_2016, 'INTERSECT', nb_lidar_2016,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2016, 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/CC_after_lidar_2016.shp')

#### 2017 ####
# CC 2017
CC_after_2017 = arcpy.SelectLayerByAttribute_management(nb_harvest_sites, 'NEW_SELECTION', '"HARVYR" > 2017')
arcpy.CopyFeatures_management(CC_after_2017, 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_after_2017.shp')
# LiDAR 2017
nb_lidar_2017 = arcpy.SelectLayerByAttribute_management(nb_tiles_index, 'NEW_SELECTION', '"Yearnr" = 2017')
arcpy.CopyFeatures_management(nb_lidar_2017, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2017.shp')

CC_after_2017 = 'Y:/William/RIP_NATURE_2/Canada/Post_harvest/CC_after_2017.shp'
nb_lidar_2017 = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/nb_lidar_2017.shp'

CC_after_lidar_2017 = arcpy.SelectLayerByLocation_management(CC_after_2017, 'INTERSECT', nb_lidar_2017,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(CC_after_lidar_2017, 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/CC_after_lidar_2017.shp')

# Merge reference sites
CC_after_lidar_2015 = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/CC_after_lidar_2015.shp'
CC_after_lidar_2016 = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/CC_after_lidar_2016.shp'
CC_after_lidar_2017 = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/CC_after_lidar_2017.shp'
# cc after 2018 had no intersections with lidar 2018
merge_nb_reference_sites = 'Y:/William/RIP_NATURE_2/Canada/dissolve/merged_buffer_within_10m_of_CC.shp'
arcpy.Merge_management([CC_after_lidar_2015, CC_after_lidar_2016, CC_after_lidar_2017], merge_nb_reference_sites)<jupyter_output><empty_output><jupyter_text>Intersect reference sites with streams<jupyter_code># inputs
merge_nb_reference_sites = 'Y:/William/RIP_NATURE_2/Canada/dissolve/merged_buffer_within_10m_of_CC.shp'
streams = 'Y:/William/RIP_NATURE_2/Canada/streams/PropertyStreams.shp'

# output
Intersected_reference_streams_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Intersected_reference_streams.shp'

# Intersect
infeatures = (streams , merge_nb_reference_sites)
arcpy.Intersect_analysis(infeatures, Intersected_reference_streams_nb)

# add geometry and select stream segments longer than 20 m.
Intersected_reference_streams_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Intersected_reference_streams.shp'
Intersected_reference_streams_over_20m_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Intersected_reference_streams_over_20m.shp'

# Add lenght
arcpy.AddGeometryAttributes_management(Intersected_reference_streams_nb, 'LENGTH', 'METERS')

# Select and copy segments over 20 m
attribute_selected_ref = arcpy.SelectLayerByAttribute_management(Intersected_reference_streams_nb, 'NEW_SELECTION', '"LENGTH" > 20')
arcpy.CopyFeatures_management(attribute_selected_ref, Intersected_reference_streams_over_20m_nb)<jupyter_output><empty_output><jupyter_text>The 20 m segments were buffered by 20 m and intersected with and dissolved by the unharvested sites.<jupyter_code>import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.parallelProcessingFactor = "100%"


Intersected_reference_streams_over_20m_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Intersected_reference_streams_over_20m.shp'
Intersected_reference_streams_over_20m_nb_Buffered = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Intersected_reference_streams_over_20m_buffered_by_20m.shp'


# Buffer stream segments by 20 m
arcpy.Buffer_analysis(Intersected_reference_streams_over_20m_nb, Intersected_reference_streams_over_20m_nb_Buffered, '20 Meters')
merge_nb_reference_sites = 'Y:/William/RIP_NATURE_2/Canada/dissolve/merged_buffer_within_10m_of_CC.shp'
Intersect_buffered_stream_segments_with_unharvested_sites_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Buffered_reference_streams_intersect_CC.shp'


# Intersect Buffered stream segments with unharvest sites
infeatures = (Intersected_reference_streams_over_20m_nb_Buffered, merge_nb_reference_sites)
arcpy.Intersect_analysis(infeatures, Intersect_buffered_stream_segments_with_unharvested_sites_nb)

# DISSOLVE BY HARVEST ID
Dissolve_reference_buffer_zones_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'
dissolveFields = ["OBJECTID_1"]
arcpy.Dissolve_management(Intersect_buffered_stream_segments_with_unharvested_sites_nb, Dissolve_reference_buffer_zones_nb, dissolveFields, "", "SINGLE_PART", "DISSOLVE_LINES")
# number of reference harvest sites

# Harvestes sites 
Harvested_nb = 'Y:/William/RIP_NATURE_2/canada/nb_buffers_lidar_post_harvest.shp'
print(arcpy.GetCount_management(Harvested_nb),'Harvested buffers')


# Reference sites
Reference_nb = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'
print(arcpy.GetCount_management(Reference_nb),'Reference unharvested buffers')


#good example location 66,2302749w 46,6845070<jupyter_output><empty_output><jupyter_text>Locate and download relevant lidar dataSelect LiDAR squares that intersects pre_harvest and post_harvest sites<jupyter_code># LiDAR tiles index
nb_tiles_index = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/geonb_li-idl_shp/geonb_li_idl_cgvd2013.shp'

# Harvested sites
Lidar_Post_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/canada/nb_buffers_lidar_post_harvest.shp'

# Reference sites
Lidar_Pre_harvest_buffer_zones = 'Y:/William/RIP_NATURE_2/Canada/Pre_harvest/Buffered_reference_streams_intersect_CC_dissolve_by_CC.shp'

# Spatial select
# Harvested
tiles_intersect_harvested = arcpy.SelectLayerByLocation_management(nb_tiles_index, 'INTERSECT', Lidar_Post_harvest_buffer_zones,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(tiles_intersect_harvested, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/geonb_li-idl_shp/tiles_intersect_harvested.shp')

# Reference
tiles_intersect_reference = arcpy.SelectLayerByLocation_management(nb_tiles_index, 'INTERSECT', Lidar_Pre_harvest_buffer_zones,'', 'NEW_SELECTION') 
arcpy.CopyFeatures_management(tiles_intersect_reference, 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/geonb_li-idl_shp/tiles_intersect_reference.shp')

## Merge
Merge_post_harvest_with_pre_harvest = 'Y:/William/RIP_NATURE_2/Canada/LiDAR_Years/geonb_li-idl_shp/Merge_post_harvest_with_pre_harvest.shp'
arcpy.Merge_management([tiles_intersect_harvested, tiles_intersect_reference], Merge_post_harvest_with_pre_harvest)<jupyter_output><empty_output><jupyter_text>Download relevant lidar files<jupyter_code>import os
from urllib.parse import urlparse
import pandas as pd
import wget
from multiprocessing.dummy import Pool 

dest_dir = 'Y:/William/RIP_NATURE_2/Canada/Download Las Files/Download Las Files/Download/'

def download_files(urls):
    wget.download(urls, dest_dir)
    
textfile = pd.read_csv('Y:/William/RIP_NATURE_2/Canada/Download Las Files/Download Las Files/Merge_post_harvest_with_pre_harvest.csv', sep =';')
remaining_files = []

url_list = textfile['FILE_URL'].tolist()
for url in url_list:
    name = os.path.basename(url)
    if name not in os.listdir(dest_dir):
        remaining_files.append(url)
print(len(remaining_files),' files are remaining.')
  
Pool(15).map(download_files,remaining_files)<jupyter_output><empty_output><jupyter_text>Extract canopy cover from lidar data lastoolslasheight -cpu64 -lof file_list.11740.txt -cores 32 -replace_z -odir "Y:\William\RIP_NATURE_2\Canada\canopy\las_height" -olazlascanopy -cpu64 -lof file_list.29684.txt -cores 32 -step 1 -height_cutoff 5 -max -odir "y:\William\RIP_NATURE_2\Canada\canopy\las_canopy" -otif Reclassify canopy so all cells over 5 m is 1 and everything else 0.<jupyter_code>import os
import sys
import time
import whitebox
wbt = whitebox.WhiteboxTools()

unclassed_canopy = 'Y:/William/RIP_NATURE_2/Canada/canopy/las_canopy/'
reclassed_canopy = 'Y:/William/RIP_NATURE_2/Canada/canopy/CanopyReclassify/'

for file in os.listdir(unclassed_canopy):
    if file.endswith('.tif') and file not in os.listdir(reclassed_canopy):
        input1 = unclassed_canopy + file
        input2 = 5
        output = reclassed_canopy + file
        wbt.greater_than(
            input1,
            input2,
            output,
            incl_equals=False
        )
print('done')

#Use copy raster to make sure no files are compressed
#Use copy raster to make sure no files are compressed
import os
import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.compression = "NONE"
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("NAD 1983 CSRS New Brunswick Stereographic")
reclassed_canopy = 'Y:/William/RIP_NATURE_2/Canada/canopy/CanopyReclassify/'
decompressed = 'Y:/William/RIP_NATURE_2/Canada/canopy/decompressed/'
for file in os.listdir(reclassed_canopy):
    if file.endswith('.tif'):
        in_raster = reclassed_canopy + file
        out_rasterdataset = decompressed + file
        arcpy.management.CopyRaster(in_raster, out_rasterdataset)<jupyter_output><empty_output><jupyter_text>Build viritual raster<jupyter_code>print('creating geodatabase for labels')
arcpy.CreateFileGDB_management('D:/GIS/RIP2/canada_canopy_mosaic/canopy_db', 'canadacanopy.gdb')
mdnameLabel = "Canada"
gdbname = "D:/GIS/RIP2/canada_canopy_mosaic/canopy_db/canadacanopy.gdb"
mdname = "canopy_mosaic"
prjfile = ""
noband = "1"
pixtype = "8_BIT_UNSIGNED"
pdef = "NONE"
wavelength = ""
print('Creating a ampty raster dataset within the geodatabase')
arcpy.CreateMosaicDataset_management(gdbname, mdname, prjfile, noband,pixtype, pdef, wavelength)

print('Adding Binary labels to dataset')
mdname = "D:/GIS/RIP2/canada_canopy_mosaic/canopy_db/canadacanopy.gdb/canopy_mosaic"
rastype = "Raster Dataset"
inpath = 'D:/GIS/RIP2/temp/reclass_canopys'
arcpy.AddRastersToMosaicDataset_management(mdname,  rastype, inpath)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0975263008245664, 'fraction_numerical': 0.06826184904092211, 'mean_word_length': 9.437507133888825, 'pattern_counts': {'":': 0, '<': 179, '<?xml version=': 0, '>': 207, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '40', 'dir': 'jupyter-structured-clean-dedup', 'id': '552402', 'license': 'MIT', 'n_tokens_mistral': 41825, 'n_tokens_neox': 36405, 'n_words': 3717, 'path': 'Salvaged_notebook_RIP2.ipynb', 'repo_name': 'williamlidberg/RIP-Nature', 'size': '91442'}