Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have been working on a Python script to automate the steps illustrated by whuber here:Create points around the edge of viewshed results

I am trying to create a horizon graph illustrating the maximum viewing angles from a single observer point in all directions. Here is the script:

import arcpy
from arcpy import env
from arcpy.sa import *

env.Workspace = "C:/Temp/Python"
env.overwriteOutput = True

inDEM = "C:/Temp/multiDEM.tif"
env.extent = inDEM

inPoints = "C:/Temp/Points_ALL.shp"
outPoints = "C:/Temp/Python/output/pts.shp"
outPath = "C:/Temp/Python/output/observers"
outputPath = "C:/Temp/Python/output"
outEuclidean = "C:/Temp/Python/output/euclidean"
fieldFID = 'FID'
Elevations = 'RASTERVALU'
arcpy.CheckOutExtension("Spatial")

ExtractValuesToPoints(inPoints, inDEM, outPoints)

arcpy.MakeFeatureLayer_management (outPoints, "pts")

with arcpy.da.SearchCursor('pts',[fieldFID, Elevations]) as cursor:
    for row in cursor:
        fid = str(row[0])
        print fid
        ### Select and export each feature
        arcpy.SelectLayerByAttribute_management ("pts", "NEW_SELECTION", '"FID" = {}'.format(fid))
        arcpy.FeatureClassToFeatureClass_conversion("pts", outPath, fid+".shp")
        ### Create Viewsheds
        outViewshed = Viewshed(inDEM,"pts",1,"CURVED_EARTH",0.15)
        outViewshed.save("C:/Temp/Python/output/viewsheds/"+fid)
        ### Euclidean Distance
        outEucDistance = EucDistance(outPath+"/"+fid+".shp")
        outEucDistance.save(outEuclidean+"/"+fid)
        Euclid = outEuclidean+"/"+fid
        ### Subtraction for Relative Elevations
        elevation = row[1]
        outMinus = Minus(inDEM, elevation)
        ### Division for "viewing angle"
        outDivide = Divide(outMinus, outEucDistance)
        ### Tangent Angles
        outTimes = Times(outViewshed, outDivide)
        outSetNull = SetNull(outTimes, outTimes, "Value = 0")
        tangentAngle = Tan(outSetNull)
        tangentAngle.save(outputPath+"/angles/"+fid)
        ### Aspect of Euclidean Distance
        outAspect = Aspect(Euclid)
        outInt = Int(outAspect)
        outInt.save("C:/Temp/Python/output/int/"+fid)
        ### Zonal Statistics
        ZonalStatisticsAsTable(outInt, "VALUE", tangentAngle, outputPath+"/tables/"+fid, "DATA", "MAXIMUM")

So this script does the following:

  1. Iterate through Points_ALL.shp to select one row at a time
  2. Export each row to a feature class
  3. Create a viewshed of each using inDEM
  4. Find the Euclidean Distance
  5. Find relative elevation through Subtraction
  6. Find Viewing angle - Dividing relative elevation by the Euclidean Distance and taking its TANgent to find the viewing angle
  7. Deriving the aspect of the Observer (then discretizing it using INT)
  8. Running ZONAL Statistics on the Integer Zones and finding the MAXIMUM angles from step 6 in each aspect "zone."

So what's the problem? Well...

After investigating each table, I noticed that I receive two errors:

  1. The maximum viewing angles in some areas are greater than 90 degrees (some as high as 600) so I think I am missing something in my calculations.
  2. Not all zones end up in the resulting table. I'll have missing aspect values even though I check the data and there are obviously results in the directions that are missing.

Can someone help me find what I am missing? Did I skip a step or misinterpret a step from the previous question?

share|improve this question
    
(1) The missing zones question has appeared before: whenever a zone contains no cell center, it won't be included in the output. In your calculation that is likely to happen, because a nearby obstruction will block a range of angles of view. (2) Before you script any geoprocessing workflow, it is a very good idea to perform it manually on several different datasets to make sure it's working correctly. (The Raster Calculator is your friend here.) Have you done that successfully? (3) You application of the tangent function to a ratio makes no sense. Perhaps you intend the inverse tangent? –  whuber May 3 '13 at 16:44
    
(1) That's what I was thinking. (2) I did this, but I think the problem is with... (3) This part. I know I misinterpreted this part. Feel free to propose a correction. I am mainly trying to interpret what you mentioned in the other question. What's the best way to fix this? –  Spatial Pariah May 3 '13 at 17:41
    
In that case, you don't have to apply the tangent at all: the ratio is the tangent of the viewing angle. –  whuber May 3 '13 at 17:58
    
Ah, that makes more sense then. Did you think that the rest looks like it "should" work? –  Spatial Pariah May 3 '13 at 20:45

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.