Skip to content

Commit f676693

Browse files
marecotecmarecotec
committed
Added Class 12
1 parent 6f6bd5a commit f676693

10 files changed

Lines changed: 340 additions & 0 deletions

File tree

‎Classes/12_Rasters/ALL_FILES.zip‎

963 KB
Binary file not shown.

‎Classes/12_Rasters/Step_0.py‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
2+
#####
3+
# Step 0 - Practice tasks before we start.
4+
#####
5+
6+
# Task a: Run the buffer tool on Step_0_Data.zip/RI_Forest_Health_Works_Project_Points_All_Invasives.shp, with a
7+
# distance of 1 mile:
8+
9+
10+
# Task b: Dissolve your resulting buffer:
11+
12+
13+
# Task c: On the original point file (RI_Forest_Health_Works_Project_Points_All_Invasives.shp), use a
14+
# search cursor to print the "Owner" field within the attributes.
15+

‎Classes/12_Rasters/Step_0_Answers.py‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
2+
#####
3+
# Step 0 - Practice tasks before we start.
4+
#####
5+
6+
# Task a: Run the buffer tool on Step_0_Data.zip/RI_Forest_Health_Works_Project_Points_All_Invasives.shp, with a
7+
# distance of 1 mile:
8+
import arcpy
9+
10+
arcpy.env.workspace = r"Z:\Andy's Documents\Teaching and students\URI\NRS - GIS Python Course\Github\Course_ArcGIS_Python\Classes\12_Rasters\Step_0_Data"
11+
12+
arcpy.Buffer_analysis("RI_Forest_Health_Works_Project_Points_All_Invasives.shp",
13+
"RI_Forest_Health_Works_Project_Points_All_Invasives_Buffer.shp",
14+
"1 mile")
15+
16+
# Task b: Dissolve your resulting buffer:
17+
arcpy.Dissolve_management("RI_Forest_Health_Works_Project_Points_All_Invasives_Buffer.shp",
18+
"RI_Forest_Health_Works_Project_Points_All_Invasives_Dissolve.shp")
19+
20+
21+
# Task c: On the original point file (RI_Forest_Health_Works_Project_Points_All_Invasives.shp), use a
22+
# search cursor to print the "Owner" field within the attributes.
23+
24+
input_shp = "RI_Forest_Health_Works_Project_Points_All_Invasives.shp"
25+
fields = ['Owner']
26+
27+
with arcpy.da.SearchCursor(input_shp, fields) as cursor:
28+
for row in cursor:
29+
print(u'Owner = {0}'.format(row[0]))
30+

‎Classes/12_Rasters/Step_0_Data.zip‎

154 KB
Binary file not shown.

‎Classes/12_Rasters/Step_1.py‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
2+
#####
3+
# Step 1 - Zonal stats and extracting values to points
4+
#####
5+
6+
# Raster processing in arcpy is very similar to working with shapefiles and other feature classes. In this example, we
7+
# conduct some zonal statistics and undertake point value extraction from a raster.
8+
9+
# Using Step_1_Data.zip, extract teh zonal values for the shapefile: Biogeography_Made_Up.shp, and the sst_mean.tif
10+
# raster.
11+
import arcpy
12+
arcpy.CheckOutExtension("Spatial")
13+
arcpy.env.overwriteOutput = True
14+
arcpy.env.workspace = r"Z:\Andy's Documents\Teaching and students\URI\NRS - GIS Python Course\Github\Course_ArcGIS_Python\Classes\12_Rasters\Step_1_Data"
15+
arcpy.gp.ZonalStatisticsAsTable_sa("Biogeography_Made_Up.shp", "Area", "sst_mean.tif", "sst_mean_zones.dbf", "DATA", "ALL")
16+
17+
# Task 1 - Using a for loop, process all three *.tif files for SST (mean, min and max), you will need to edit the code above.
18+
# HInt: ZonalStats... is sensitive to file name output, you may need to use something that shortebs the output table name:
19+
# e.g. file_name = raster[:6] + "_out.dbf"
20+
21+
22+
23+
# Pulling values from a raster using a point shapefule is also pretty easy, this time, we use the Extract Values tool.
24+
25+
arcpy.gp.ExtractValuesToPoints_sa("Great_Whites.shp", "sst_mean.tif", "Great_Whites_Extract.shp", "NONE", "VALUE_ONLY")
26+
27+
# Task 2 - I want you to run the tool above on the rasters provided (mean, min and max). Think about how you will do this,
28+
# maybe search online about extracting multiple values to points.
29+

‎Classes/12_Rasters/Step_1_Answers.py‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
2+
#####
3+
# Step 1 - Zonal stats and extracting values to points
4+
#####
5+
6+
# Raster processing in arcpy is very similar to working with shapefiles and other feature classes. In this example, we
7+
# conduct some zonal statistics and undertake point value extraction from a raster.
8+
9+
# Using Step_1_Data.zip, extract teh zonal values for the shapefile: Biogeography_Made_Up.shp, and the sst_mean.tif
10+
# raster.
11+
import arcpy
12+
arcpy.CheckOutExtension("Spatial")
13+
arcpy.env.overwriteOutput = True
14+
arcpy.env.workspace = r"Z:\Andy's Documents\Teaching and students\URI\NRS - GIS Python Course\Github\Course_ArcGIS_Python\Classes\12_Rasters\Step_1_Data"
15+
arcpy.gp.ZonalStatisticsAsTable_sa("Biogeography_Made_Up.shp", "Area", "sst_mean.tif", "sst_mean_zones.dbf", "DATA", "ALL")
16+
17+
# Task 1 - Using a for loop, process all three *.tif files for SST (mean, min and max), you will need to edit the code above.
18+
# HInt: ZonalStats... is sensitive to file name output, you may need to use something that shortebs the output table name:
19+
# e.g. file_name = raster[:6] + "_out.dbf"
20+
21+
raster_list = arcpy.ListRasters()
22+
23+
for raster in raster_list:
24+
file_name = raster[:6] + "_out.dbf"
25+
arcpy.gp.ZonalStatisticsAsTable_sa("Biogeography_Made_Up.shp", "Area", raster, file_name, "DATA",
26+
"ALL")
27+
28+
# Pulling values from a raster using a point shapefule is also pretty easy, this time, we use the Extract Values tool.
29+
30+
arcpy.gp.ExtractValuesToPoints_sa("Great_Whites.shp", "sst_mean.tif", "Great_Whites_Extract.shp", "NONE", "VALUE_ONLY")
31+
32+
# Task 2 - I want you to run the tool above on the rasters provided (mean, min and max). Think about how you will do this,
33+
# maybe search online about extracting multiple values to points.
34+
35+
raster_list = arcpy.ListRasters()
36+
arcpy.gp.ExtractMultiValuesToPoints("Great_Whites.shp", raster_list, "BILINEAR")

‎Classes/12_Rasters/Step_1_Data.zip‎

778 KB
Binary file not shown.

‎Classes/12_Rasters/Step_2.py‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
2+
#####
3+
# Step 2 - Interpolation and routines to conduct validation
4+
#####
5+
6+
# We will interpolate a file of elevations, and use a k-fold validation approach, whereby we will loop through creating
7+
# n interpolations, validate with points that were dropped from the interpolation, and report
8+
# a correlation value using Pandas.
9+
10+
# This is a multi-step process:
11+
12+
# 1) Let's do a random selection, and save the split files in different directories.
13+
14+
import arcpy
15+
import random as random
16+
from arcpy.sa import *
17+
arcpy.CheckOutExtension("Spatial")
18+
from scipy.stats.stats import pearsonr
19+
20+
arcpy.env.workspace = r"Z:\Andy's Documents\Teaching and students\URI\NRS - GIS Python Course\Github\Course_ArcGIS_Python\Classes\12_Rasters\Step_2_Data"
21+
input_data = r"Elevations.shp"
22+
23+
arcpy.env.extent = input_data
24+
25+
kfold = 20 #number of folds
26+
prop_samples = 30 #proportion of points to split
27+
i = 0
28+
29+
30+
# Code Obtained from: https://support.esri.com/en/technical-article/000013141, edited to return 2 outputs
31+
def SelectRandomByPercent (layer, layer2, percent):
32+
fc = arcpy.Describe(layer).catalogPath
33+
featureCount = float(arcpy.GetCount_management (fc).getOutput(0))
34+
count = int(featureCount * float(percent) / float(100))
35+
if not count:
36+
arcpy.SelectLayerByAttribute_management (layer, "CLEAR_SELECTION")
37+
return
38+
oids = [oid for oid, in arcpy.da.SearchCursor(fc, "OID@")]
39+
oidFldName = arcpy.Describe(layer).OIDFieldName
40+
delimOidFld = arcpy.AddFieldDelimiters (layer, oidFldName)
41+
randOids = random.sample (oids, count)
42+
oidsStr = ", ".join(map(str, randOids))
43+
sql = "{0} IN ({1})".format(delimOidFld, oidsStr)
44+
output1 = arcpy.SelectLayerByAttribute_management(layer, "", sql)
45+
sql = "{0} NOT IN ({1})".format(delimOidFld, oidsStr)
46+
output2 = arcpy.SelectLayerByAttribute_management(layer2, "", sql)
47+
return output1, output2
48+
49+
50+
while i < kfold:
51+
arcpy.MakeFeatureLayer_management(input_data, "input_data_lyr")
52+
arcpy.MakeFeatureLayer_management(input_data, "input_data_lyr2")
53+
output1, output2 = SelectRandomByPercent("input_data_lyr", "input_data_lyr2", prop_samples)
54+
arcpy.CopyFeatures_management(output1, "elv_" + str(i) + "_test.shp")
55+
arcpy.CopyFeatures_management(output2, "elv_" + str(i) + "_train.shp")
56+
arcpy.Delete_management("input_data_lyr")
57+
arcpy.Delete_management("input_data_lyr2")
58+
i = i + 1
59+
60+
61+
# 2) Now let's create several surfaces using simple interpolation:
62+
63+
training_shp = arcpy.ListFeatureClasses("*train*")
64+
count = 0
65+
66+
for i in training_shp:
67+
outIDW = Idw(i, "Elevation", 2000, 2, RadiusVariable(10, 50000))
68+
outIDW.save("elev_" + str(count) + ".tif")
69+
count = count + 1
70+
71+
72+
73+
# Task - Now we have run through the generation and partitioning of the input shapefile, and created a
74+
# rudimentary interpolation using IDW. We need to validate the file. You should consider the steps you need to
75+
# undertake to validate each partition. I have given you a run down below:
76+
77+
# 1. You need to list the "test" files, see line 58...
78+
# 2. For each IDW and each test file (they share numbers... so you can use a count variable to iterate this,
79+
# you need to use ExtractValuesToPoints (in_point_features, in_raster, out_point_features) to pull the test
80+
# values from the "test" shapefile.
81+
# 3. Use the code I provide below to get yourself a correlation value for the interpolation.
82+
83+
# arr = arcpy.da.FeatureClassToNumPyArray("elev_" + str(count) + "_valid.shp", ["Elevation", "RASTERVALU"])
84+
# print "Interpolation " + str(count) + " correlation = " + str(pearsonr(arr["Elevation"], arr["RASTERVALU"]))
85+
86+
# 4. Store the correlation value (hint the pearsonr outputs a tuple), take the mean of the total k-fold (hint sum / length).
87+
88+
# 5. Take a mean of all the interpolated surfaces you generated and store it as an output (Hint:
89+
# https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/cell-statistics.htm
90+
91+
# 6. Take a standard deviation of all the interpolated surfaces you generated and store it as an output (Hint:
92+
# https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/cell-statistics.htm
93+
94+
# 6. Add clean up code to remove temporary files.
95+

‎Classes/12_Rasters/Step_2_Answers.py‎

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
2+
#####
3+
# Step 2 - Interpolation and routines to conduct validation
4+
#####
5+
6+
# We will interpolate a file of elevations, and use a k-fold validation approach, whereby we will loop through creating
7+
# n interpolations, validate with points that were dropped from the interpolation, and report
8+
# a correlation value using Pandas.
9+
10+
# This is a multi-step process:
11+
12+
# 1) Let's do a random selection, and save the split files in different directories.
13+
14+
import arcpy
15+
import random as random
16+
from arcpy.sa import *
17+
arcpy.CheckOutExtension("Spatial")
18+
from scipy.stats.stats import pearsonr
19+
20+
arcpy.env.workspace = r"Z:\Andy's Documents\Teaching and students\URI\NRS - GIS Python Course\Github\Course_ArcGIS_Python\Classes\12_Rasters\Step_2_Data"
21+
input_data = r"Elevations.shp"
22+
23+
arcpy.env.extent = input_data
24+
25+
kfold = 20 #number of folds
26+
prop_samples = 30 #proportion of points to split
27+
i = 0
28+
29+
30+
# Code Obtained from: https://support.esri.com/en/technical-article/000013141, edited to return 2 outputs
31+
def SelectRandomByPercent (layer, layer2, percent):
32+
fc = arcpy.Describe(layer).catalogPath
33+
featureCount = float(arcpy.GetCount_management (fc).getOutput(0))
34+
count = int(featureCount * float(percent) / float(100))
35+
if not count:
36+
arcpy.SelectLayerByAttribute_management (layer, "CLEAR_SELECTION")
37+
return
38+
oids = [oid for oid, in arcpy.da.SearchCursor(fc, "OID@")]
39+
oidFldName = arcpy.Describe(layer).OIDFieldName
40+
delimOidFld = arcpy.AddFieldDelimiters (layer, oidFldName)
41+
randOids = random.sample (oids, count)
42+
oidsStr = ", ".join(map(str, randOids))
43+
sql = "{0} IN ({1})".format(delimOidFld, oidsStr)
44+
output1 = arcpy.SelectLayerByAttribute_management(layer, "", sql)
45+
sql = "{0} NOT IN ({1})".format(delimOidFld, oidsStr)
46+
output2 = arcpy.SelectLayerByAttribute_management(layer2, "", sql)
47+
return output1, output2
48+
49+
50+
while i < kfold:
51+
arcpy.MakeFeatureLayer_management(input_data, "input_data_lyr")
52+
arcpy.MakeFeatureLayer_management(input_data, "input_data_lyr2")
53+
output1, output2 = SelectRandomByPercent("input_data_lyr", "input_data_lyr2", prop_samples)
54+
arcpy.CopyFeatures_management(output1, "elv_" + str(i) + "_test.shp")
55+
arcpy.CopyFeatures_management(output2, "elv_" + str(i) + "_train.shp")
56+
arcpy.Delete_management("input_data_lyr")
57+
arcpy.Delete_management("input_data_lyr2")
58+
i = i + 1
59+
60+
61+
# 2) Now let's create several surfaces using simple interpolation:
62+
63+
training_shp = arcpy.ListFeatureClasses("*train*")
64+
count = 0
65+
66+
for i in training_shp:
67+
outIDW = Idw(i, "Elevation", 2000, 2, RadiusVariable(10, 50000))
68+
outIDW.save("elev_" + str(count) + ".tif")
69+
count = count + 1
70+
71+
72+
73+
# Task - Now we have run through the generation and partitioning of the input shapefile, and created a
74+
# rudimentary interpolation using IDW. We need to validate the file. You should consider the steps you need to
75+
# undertake to validate each partition. I have given you a run down below:
76+
77+
# 1. You need to list the "test" files, see line 58...
78+
# 2. For each IDW and each test file (they share numbers... so you can use a count variable to iterate this,
79+
# you need to use ExtractValuesToPoints (in_point_features, in_raster, out_point_features) to pull the test
80+
# values from the "test" shapefile.
81+
# 3. Use the code I provide below to get yourself a correlation value for the interpolation.
82+
83+
# arr = arcpy.da.FeatureClassToNumPyArray("elev_" + str(count) + "_valid.shp", ["Elevation", "RASTERVALU"])
84+
# print "Interpolation " + str(count) + " correlation = " + str(pearsonr(arr["Elevation"], arr["RASTERVALU"]))
85+
86+
# 4. Store the correlation value (hint the pearsonr outputs a tuple), take the mean of the total k-fold (hint sum / length).
87+
88+
# 5. Take a mean of all the interpolated surfaces you generated and store it as an output (Hint:
89+
# https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/cell-statistics.htm
90+
91+
# 6. Take a standard deviation of all the interpolated surfaces you generated and store it as an output (Hint:
92+
# https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/cell-statistics.htm
93+
94+
# 6. Add clean up code to remove temporary files.
95+
96+
test_shp = arcpy.ListFeatureClasses("*test*")
97+
pearson_values = []
98+
raster_list = []
99+
100+
for fc in test_shp:
101+
run = fc.split("_")
102+
arcpy.gp.ExtractValuesToPoints_sa(fc, "elev_" + str(run[1]) + ".tif", "elev_" + str(run[1]) + "_valid.shp", "NONE", "VALUE_ONLY")
103+
arr = arcpy.da.FeatureClassToNumPyArray("elev_" + str(run[1]) + "_valid.shp", ["Elevation", "RASTERVALU"])
104+
print "Interpolation " + str(run[1]) + " correlation = " + str(pearsonr(arr["Elevation"], arr["RASTERVALU"]))
105+
106+
pearson_values.append(pearsonr(arr["Elevation"], arr["RASTERVALU"])[0])
107+
raster_list.append("elev_" + str(run[1]) + ".tif")
108+
109+
arcpy.Delete_management("elev_" + str(run[1]) + "_valid.shp")
110+
arcpy.Delete_management("elv_" + str(run[1]) + "_train.shp")
111+
arcpy.Delete_management("elv_" + str(run[1]) + "_test.shp")
112+
113+
print "Overall mean: " + str(sum(pearson_values) / float(len(pearson_values)))
114+
115+
outCellStats = arcpy.sa.CellStatistics(raster_list, "MEAN", "DATA")
116+
outCellStats.save("out_surface_mean.tif")
117+
118+
outCellStats = arcpy.sa.CellStatistics(raster_list, "STD", "DATA")
119+
outCellStats.save("out_surface_std.tif")
120+
121+
for i in raster_list:
122+
arcpy.Delete_management(i)
123+
124+
125+
126+
127+
128+
129+
130+
131+
132+
133+
134+
135+

‎Classes/12_Rasters/Step_2_Data.zip‎

25 KB
Binary file not shown.

0 commit comments

Comments
 (0)