-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGIS-featuresExtract_Trials.txt
76 lines (57 loc) · 2.96 KB
/
GIS-featuresExtract_Trials.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
fields = ['OBJECTID']
with arcpy.da.SearchCursor(fields) as cursor:
for row in cursor:
select = 'OBJECTID = {}'.format(row[0])
arcpy.SelectLayerByAttribute_management('tractsLyr','NEW_SELECTION',select)
###########################################L A T E S T############################################
##################################################################################################
#Making watershed for every point
from arcpy.sa import *
PointLayer="Points"
fields = ['OID']
#Limting number of rows in the field
exp='OID<10'
# Set local variables
inFlowDirection = "FlowDir_dem_2"
with arcpy.da.SearchCursor(PointLayer, fields, exp) as cursor:
for row in cursor:
##Select a row of PointLayer
inPourPointData = 'OID = {}'.format(row[0])
arcpy.SelectLayerByAttribute_management(PointLayer,'NEW_SELECTION',select)
##Create a temporary feature class
####Temporary point layer
out_path = "D:/VU/Re/DC_Data"
out_name = "temp-pour.shp"
geometry_type = "POINT"
template = "Points"
has_m = "DISABLED"
has_z = "DISABLED"
# Use Describe to get a SpatialReference object
spatial_ref = arcpy.Describe("D:/VU/Re/DC_Data/dem_ratertopoint.shp").spatialReference
# Execute CreateFeatureclass
arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template,
has_m, has_z, spatial_ref)
###Copy from original point layer
arcpy.management.CopyFeatures(inPourPointData, out_name)
###Execute Watershed
outWatershed = Watershed('FlowDir_Fill1', out_name)
# Save the output
inR = outWatershed.save("D:/VU/Re/DC_Data/outwtrshd02.tif")
###We need watershed polygon layer...
outpoly = "D:/VU/Re/DC_Data/outwtr.shp"
arcpy.RasterToPolygon_conversion(outWatershed, outpoly)
arcpy.SelectLayerByAttribute_management(outpoly,'NEW_SELECTION')
# Set local variables
inPointFeatures = arcpy.SelectLayerByLocation_management(PointLayer,"INTERSECT",outpoly,"","NEW_SELECTION")
inRaster = "Slope_Fill_d3"
outPointFeatures = "C:/Users/Himel/Documents/ArcGIS/Projects/MyProject/extractvaluespts.shp"
# Execute ExtractValuesToPoints ###Why not extracting only the rastervalue field???
ExtractValuesToPoints(inPointFeatures, inRaster, outPointFeatures,
"INTERPOLATE", "VALUE_ONLY")
###Step 3: Calculate the average and add in the main watershed table
na = arcpy.da.TableToNumPyArray('extractvaluespts', 'RASTERVALU')
mean = numpy.mean(na['RASTERVALU'])
arcpy.management.DeleteFeatures(out_name)
print(mean)
####################################################################################
####################################################################################