POST
|
Hi Xander, Thanks for your help so far and sorry for my ignorance. I tried converting my points to feature classes and used a fgdb raster, then tried your initial code again, but got the same errors. The points overlap with the raster and are in the same coordinate system, and I set Union of Inputs as the processing extent (though do I actually need to do this in Python?). The vast majority of cells have values. I have attached the TIFF file from which the fgdb raster was made. >>> import arcpy, os, numpy
... from numpy import *
...
... # settings of your data
... fcpoints = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clip16points'
... fldVisible = 'VisiblePix' # will contain the number of visible pixels
... dem = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clipdem30mWGS1'
...
... # ouput workspace and name for viewshed
... outputWS = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb'
... prefixVS = 'Viewshed'
... arcpy.env.workspace = outputWS
... arcpy.env.overwriteOutput = True
...
... # grab a SA license
... if arcpy.CheckExtension("Spatial") == "Available":
... arcpy.CheckOutExtension("Spatial")
... else:
... exit("no SA license available")
... desc = arcpy.Describe(fcpoints)
... if desc.hasOID:
... fldOID = desc.OIDFieldName
... else:
... exit("featureclass does not have an Object ID field")
...
... # loop through points
... fields = ("OID@",fldVisible)
... with arcpy.da.UpdateCursor(fcpoints, fields) as cursor:
... for row in cursor:
... oid = row[0]
...
... # define output name (of oid = 1, raster will be called Viewshed1
... outVSname = outputWS + os.sep + prefixVS + str(oid)
...
... # create featurelayer for current observer point
... whereclause = "fldOID = {0}".format(oid)
... FLsources = "sources_Layer"
... arcpy.MakeFeatureLayer_management(fcpoints,FLsources,whereclause,outputWS)
...
... # calculate viewshed
... arcpy.gp.Viewshed_sa(dem,FLsources,outVSname,"1","FLAT_EARTH","0,13")
...
... # determine number of visible cells
... myArray = arcpy.RasterToNumPyArray(outVSname)
... myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape))
... visPixels = myArray2.sum()
... print "There are {0} visible pixels for point with OID={1}".format(visPixels,oid)
...
... # update value in point featureclass
... row[1] = visPixels
... cursor.updateRow(row)
... del myArray, myArray2
...
... del row
...
... # release license
... arcpy.CheckInExtension("Spatial") The processing output is also attached. [HTML]Executing: MakeFeatureLayer C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clip16points sources_Layer "fldOID = 1" C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb "OBJECTID_1 OBJECTID VISIBLE NONE;Shape Shape VISIBLE NONE;Join_Count Join_Count VISIBLE NONE;TARGET_FID TARGET_FID VISIBLE NONE;FID_ FID_ VISIBLE NONE;Date Date VISIBLE NONE;Date_2 Date_2 VISIBLE NONE;Year Year VISIBLE NONE;Field_comp Field_comp VISIBLE NONE;Unique_ID Unique_ID VISIBLE NONE;X X VISIBLE NONE;Y Y VISIBLE NONE;Waypoint_I Waypoint_I VISIBLE NONE;elevation elevation VISIBLE NONE;Observers Observers VISIBLE NONE;Photo_ID Photo_ID VISIBLE NONE;Sward_heig Sward_heig VISIBLE NONE;Sward_he_1 Sward_he_1 VISIBLE NONE;Sward_he_2 Sward_he_2 VISIBLE NONE;Crop_type Crop_type VISIBLE NONE;General_co General_co VISIBLE NONE;Notes Notes VISIBLE NONE;Turbine_me Turbine_me VISIBLE NONE;Anser_drop Anser_drop VISIBLE NONE;RbG_droppi RbG_droppi VISIBLE NONE;Unknown_dr Unknown_dr VISIBLE NONE;Total_drop Total_drop VISIBLE NONE;Height_at_ Height_at_ VISIBLE NONE;Height_at1 Height_at1 VISIBLE NONE;Height_a_1 Height_a_1 VISIBLE NONE;Height_a_2 Height_a_2 VISIBLE NONE;Height_a_3 Height_a_3 VISIBLE NONE;Height_a_4 Height_a_4 VISIBLE NONE;Height_a_5 Height_a_5 VISIBLE NONE;Height_a_6 Height_a_6 VISIBLE NONE;Height_a_7 Height_a_7 VISIBLE NONE;Height_a_8 Height_a_8 VISIBLE NONE;Average_he Average_he VISIBLE NONE;Cat_No__Pl Cat_No__Pl VISIBLE NONE;Cat_No___1 Cat_No___1 VISIBLE NONE;Cat_No___2 Cat_No___2 VISIBLE NONE;Cat_No___3 Cat_No___3 VISIBLE NONE;Cat_No___4 Cat_No___4 VISIBLE NONE;Cat_No___5 Cat_No___5 VISIBLE NONE;Cat_No___6 Cat_No___6 VISIBLE NONE;Cat_No___7 Cat_No___7 VISIBLE NONE;Cat_No___8 Cat_No___8 VISIBLE NONE;Cat_No___9 Cat_No___9 VISIBLE NONE;No__densit No__densit VISIBLE NONE;No__dens_1 No__dens_1 VISIBLE NONE;No__dens_2 No__dens_2 VISIBLE NONE;No__dens_3 No__dens_3 VISIBLE NONE;No__dens_4 No__dens_4 VISIBLE NONE;Crop_colou Crop_colou VISIBLE NONE;Crop_col_1 Crop_col_1 VISIBLE NONE;Crop_col_2 Crop_col_2 VISIBLE NONE;Snow_Thin Snow_Thin VISIBLE NONE;Snow_Mediu Snow_Mediu VISIBLE NONE;Snow_Thick Snow_Thick VISIBLE NONE;_snow_cov _snow_cov VISIBLE NONE;Cat_dist_n Cat_dist_n VISIBLE NONE;Nearest_wa Nearest_wa VISIBLE NONE;Nearest__1 Nearest__1 VISIBLE NONE;Nearest__2 Nearest__2 VISIBLE NONE;Nearest__3 Nearest__3 VISIBLE NONE;Nearest__4 Nearest__4 VISIBLE NONE;No_visible No_visible VISIBLE NONE;Dist_trees Dist_trees VISIBLE NONE;Dist_turbi Dist_turbi VISIBLE NONE;Dist_roads Dist_roads VISIBLE NONE;WF1km WF1km VISIBLE NONE;WF500m WF500m VISIBLE NONE;WF250m WF250m VISIBLE NONE;WF2km WF2km VISIBLE NONE;Dist_track Dist_track VISIBLE NONE;Dist_MnRd Dist_MnRd VISIBLE NONE;Trees_1km Trees_1km VISIBLE NONE;km_roads km_roads VISIBLE NONE;km_trees km_trees VISIBLE NONE;m_track m_track VISIBLE NONE;m_roads m_roads VISIBLE NONE;km_roads_1 km_roads_1 VISIBLE NONE;m_track_1 m_track_1 VISIBLE NONE;m_roads_1 m_roads_1 VISIBLE NONE;m_trees m_trees VISIBLE NONE;m_trees_1 m_trees_1 VISIBLE NONE;km_trees_1 km_trees_1 VISIBLE NONE;km_tracks km_tracks VISIBLE NONE;km_tracks_ km_tracks_ VISIBLE NONE;Offset_A Offset_A VISIBLE NONE;VisiblePix VisiblePix VISIBLE NONE" Start Time: Tue Dec 03 12:49:04 2013 Succeeded at Tue Dec 03 12:49:04 2013 (Elapsed Time: 0.00 seconds) Executing: Viewshed C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clipdem30mWGS1 sources_Layer C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\Viewshed1 1 FLAT_EARTH 0.13 Start Time: Tue Dec 03 12:49:08 2013 ERROR 999999: Error executing function. Cannot acquire a lock. Cannot acquire a lock. Cannot acquire a lock. Cannot acquire a lock. ERROR 010004: All cells in Raster C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\clipdem30mWGS1 have the NODATA value. Stop execution. No valid Observation points present. ERROR 010160: Unable to open raster C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\Viewshed1. ERROR 010160: Unable to open raster C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\Viewshed1. Failed to execute (Viewshed). Failed at Tue Dec 03 12:49:08 2013 (Elapsed Time: 0.00 seconds)[/HTML] Any ideas? Anne
... View more
12-03-2013
02:53 AM
|
0
|
0
|
371
|
POST
|
Thanks for your reply. I am rather new to python but have tried to apply your script to my data. However, when I run it I get the error messages: ERROR 999999: Error executing function. Cannot acquire a lock. Cannot acquire a lock. Cannot acquire a lock. Cannot acquire a lock. ERROR 010004: All cells in Raster G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\combraster2 have the NODATA value. Stop execution. No valid Observation points present. ERROR 010160: Unable to open raster C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\Viewshed0. ERROR 010160: Unable to open raster C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb\Viewshed0. Failed to execute (Viewshed). Don't understand why it says there is no data in the raster and no valid obs points present. Any ideas? My script is: >>> import arcpy, os, numpy ... from numpy import* ... ... # settings of your data ... fcpoints = r'G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\clip16points.shp' ... fldVisible = 'VisiblePix' # will contain the number of visible pixels ... dem = r'G:\xxx\xxxx\xxxx\xxxx\xxx\xxxx\xxx\combraster2' ... ... # ouput workspace and name for viewshed ... outputWS = r'C:\Users\Anne.Harrison\Documents\ArcGIS\Default2.gdb' ... prefixVS = 'Viewshed' ... arcpy.env.workspace = outputWS ... arcpy.env.overwriteOutput = True ... ... # grab a SA license ... if arcpy.CheckExtension("Spatial") == "Available": ... arcpy.CheckOutExtension("Spatial") ... else: ... exit("no SA license available") ... ... # loop through points ... fields = ("oid@",fldVisible) ... with arcpy.da.UpdateCursor(fcpoints, fields) as cursor: ... for row in cursor: ... oid = row[0] ... ... # define output name (of oid = 1, raster will be called Viewshed1 ... outVSname = outputWS + os.sep + prefixVS + str(oid) ... ... # create featurelayer for current observer point ... whereclause = "OBJECTID = {0}".format(oid) ... FLsources = "sources_Layer" ... arcpy.MakeFeatureLayer_management(fcpoints,FLsources,whereclause,outputWS) ... ... # calculate viewshed ... arcpy.gp.Viewshed_sa(dem,FLsources,outVSname,"1","FLAT_EARTH","0.13") ... ... # determine number of visible cells ... myArray = arcpy.RasterToNumPyArray(outVSname) ... myArray2 = where(myArray==1,ones(myArray.shape),zeros(myArray.shape)) ... VisiblePix = myArray2.sum() ... print "There are {0} visible pixels for point with oid={1}".format(VisiblePix,oid) ... ... # update value in point featureclass ... row[1] = VisiblePix ... cursor.updateRow(row) ... del myArray, myArray2 ... ... del row ... ... # release license ... arcpy.CheckInExtension("Spatial")
... View more
11-29-2013
06:35 AM
|
0
|
0
|
371
|
POST
|
Hi, I have used iteration in ModelBuilder to create multiple viewsheds for each observer location in a point feature class. However, I now need to transfer the value for visible pixels back to the observer feature class, so that I can then compare the total visible areas of different observer points. I can't join the tables as there is no common field in the raster and points files. Does anyone know a way of automating this without manually having to insert the values into a new field in the points attributes table? Thanks!
... View more
11-25-2013
07:03 AM
|
0
|
6
|
2090
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:23 AM
|