2

I'd like to do a viewshed analysis of mines in the midwest. What I have are point features indicating the location of the mines. I also have an acreage field, which I could use to expand from the point and create a buffer polygon that approximates the size of the mine. I'd like to do a viewshed analysis on this polygon, rather than the point, since the point is just an arbitrary "center" of the mine.

Unfortunately, Viewshed in ArcGIS only allows me to approximate the viewshed of a point itself. Are there any ways to do polygons instead?

  • You could use the vertices of the polygons as points and then mosaic the resulting viewsheds. – Tom Dilts Apr 20 '16 at 18:19
  • Thanks, everyone! This was great. I'm running it now after using Feature Vertices to Points tool. Since I'm interested in understanding views of each mine (and not all together), I'm using an iterative feature selection approach based on the original ID. – Andy B May 5 '16 at 19:21
2

If you were to convert your buffers to points (points along the perimeter) using the Feature Vertices to points tool, create the viewshed raster for each point, then use raster math to add all the rasters together, you would have a raster where the highest value is "visible", lower values are "partially visible", and 0 would be "not visible". Or you could reclassify the raster where the highest value is 1 and everything else is 0, so you would have a viewshed that would show which areas are completely visible from the polygon.

If you don't have an advanced license to convert feature vertices to points, here's an example of how to convert vertices to points using the shape object in Python. This will convert the centroid and all vertices in the polygon to point geometries in a new point feature class.

import arcpy
import json

workspace = r'your\workspace\here.gdb'
arcpy.env.workspace = workspace

polygons = "buffers"
polygonFields = [u'OBJECTID', u'SHAPE@']

spatialRef = arcpy.Describe(polygons).spatialReference

arcpy.CreateFeatureclass_management(workspace, polygons + "_points", "POINT", '','','', spatialRef)

points = polygons + "_points"
arcpy.AddField_management(points, "Buffer_ID", "LONG")
pointFields = [u'Buffer_ID', u'SHAPE@']

insert = arcpy.da.InsertCursor(points, pointFields)
print "Creating vertices to points for {}".format(polygons)

total = int(arcpy.GetCount_management(polygons).getOutput(0))
with arcpy.da.SearchCursor(polygons, polygonFields) as cursor:
    for row in cursor:
        [x,y] = row[1].centroid
        ptGeom = arcpy.PointGeometry(arcpy.Point(x, y))
        polyID = row[0]
        insert.insertRow((polyID, ptGeom))
        ringlist = json.loads(row[1].JSON)["rings"]
        for ringIndex, ring in enumerate(ringlist):
            vertexList = json.loads(row[1].JSON)["rings"][ringIndex]
            for index, vertex in enumerate(vertexList):
                [x,y] = vertex
                ptGeom = arcpy.PointGeometry(arcpy.Point(x, y))
                polyID = row[0]
                insert.insertRow((polyID, ptGeom))
del cursor

EDIT: If you multiply the rasters you wouldn't have to perform the additional step to reclassify (assuming the "partially visible" information is not useful).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.