Visibility tools for ArcGIS Desktop?

Question

I've been trying to create visibility areas for CCTV in my city. Obviously simple geoprocessing tools aren't the answer, cause they create results as shown on the photo. CCTV can't look through buildings, so I was looking for some solution that will skip areas beyond building. Is there any tool or extension in ArcGIS Desktop?

It would be also very helpful to be able to create oriented buffers (since not all cameras work in 360degrees), but I haven't found any working solution yet.

Answer 1

Building onto the solution that @GBG gave, then viewshed tool would be a good starting point for your application. Yes, the viewshed tool does need a DEM, but in this case, as everything is relative, the source of the DEM is inconsequential. Please see my results and explanation below:

Here is my result, using the viewshed tool and a mock dataset:

The Green dots are potential camera positions, the green areas are the visible areas.

I set a buffer distance of 100m and for one of the cameras I limited the view for only areas from the NE to NW (which is the middle camera in the image below - a mistake). From this output, I can convert it to Polygon (vector) and then process it as I need.

Here are the steps.

1. First create a polygon that will be your study area, a simple rectangle encompassing where you want to erun your analysis. Create a field in this layer called Value, assign a value of 1 to this field. Convert this layer to raster, using the value field.
2. Next, get your building data (obstructions) and create a field called value in this feature. Here you can assign a value larger than 1. I would suggest (if you want to be really accurate, you assign the building height relative to the ground level). It doesn't matter really, but a camera may "see" over lower buildings. For testing I randomly chose 5m. Convert this layer to raster, using the building height field.

You will now have two rasters.

3. Use the "Mosaic to new raster" Tool and combine your two rasters into one. You need to make sure that the "Mosaic Operator" is set to the position of the building height raster. So, if the building raster is FIRST set the mosaic operator to FIRST. If it is LAST in the list, then set the operator to last. Run the tool.
4. Once you have this raster, use the cross section tool to confirm that your raster is showing the elevations correctly. See my example below:

5. Next take place the positions of the cameras. If they are existing point files that is fine. You need to set the Camera layer up correctly to run the viewshed tool. In your Camera layer you need to add the following fields and values.

These fields are what you need:

• OFFSETA: This field governs the height of the camera, use the elevation above ground level. You can add this height per camera, in case each camera is set up differently.
• AZIMUTH1: This is the starting point for determining the viewing angle of the camera. Think of it as an arc, this is the starting point of the arc, the value is in degrees where 0 is north and 90 is East.
• AZIMUTH2: This is where the arc will end.

If the camera can see 360 degrees, Azimuth1 will be 0 and Azimuth2 will be 360.

• Last field is RADIUS2: this is the maximum viewing distance.

Now run the Viewshed tool, use the mosaicked raster as the input and the cameras as the input points.

You should get a result you can use.

Answer 2

The Spatial Analyst tools in ArcToolbox include a function called Viewshed.

You should be able to attribute your camera locations with camera heights and angles and get back the visible areas assuming you have a sufficient surface model.

Answer 3

Vector approach using standard (advanced license) tools:

arcpy.Buffer_analysis("centres", "../BUFFER.shp", "360")
arcpy.Densify_edit("BUFFER", "DISTANCE", "6.28333333333333 Meters")
arcpy.FeatureVerticesToPoints_management("BUFFER", "../all_points.shp", "ALL")
arcpy.CalculateField_management("all_points", "ORIG_FID", "[FID]")
arcpy.AddGeometryAttributes_management("centres", "POINT_X_Y_Z_M")
arcpy.SpatialJoin_analysis("all_points", "centres", "../manyCenters.shp", "JOIN_ONE_TO_ONE", "KEEP_ALL", match_option="CLOSEST")
arcpy.CalculateField_management(in_table="manyCenters", field="Shape", expression="arcpy.Point( !POINT_X!, !POINT_Y!)", "PYTHON_9.3")
arcpy.Merge_management("manyCenters;all_points", "../merged.shp")
arcpy.PointsToLine_management("merged", "../longLines.shp", "ORIG_FID")
arcpy.Erase_analysis("longLines", "HOUSES", "../erased.shp")
arcpy.MultipartToSinglepart_management("erased", "../s_parts.shp")
arcpy.SelectLayerByLocation_management("s_parts", "INTERSECT", "centres", selection_type="NEW_SELECTION", invert_spatial_relationship="INVERT")
arcpy.DeleteFeatures_management("s_parts")
arcpy.FeatureVerticesToPoints_management("s_parts", "../ends.shp", "END")
arcpy.PointsToLine_management("ends", "../outline.shp")

OUTPUT:

Note: in order to improve results, reduce step in Densify tool.

Answer 4

If you are really fixed on the 2-D approach you could generate lines that radiate from the camera location. The script below will take an existing line shapefile as an input, the number of lines required (you could create 360 lines if you wanted), and the coordinates of the camera as inputs. After the lines are created erase the building footprints from the lines. Ensure you have single part geometry, then select only those lines that intersect your camera. This will give you a rough idea of what is visible from the camera in 2-D space. Here is an example of what the script will output with the parameters set as shown in the script.

import sys, traceback, math, arcpy
def PolarToCartesian(r, theta):
    import math
    x = r * math.cos(theta)
    y = r * math.sin(theta)
    return [x, y]

####    User Parameters
#An existing shapefile in a projected CRS...
inFC = r"Z:\GISpublic\GerryG\PythonScripts\RadialLines\RadialLines_Project.shp"
#The location of the point of origin in a projected CRS...
orginX = -14996395
orginY = 5609063
#The number of lines to create...
lineCount = 12
#The distance of the line from the origin point...
lineDistance = 10000000.0
#The number of degrees to rotate the output lines counter clockwise from grid east...
lineShiftInDegrees = 0
#### End User Parameters


angleBetweenLines = 360.0/lineCount
radiansBetweenLines = angleBetweenLines *math.pi/180.0
lineRadians = 0.0 + lineShiftInDegrees*math.pi/180.0
cursor = arcpy.da.InsertCursor(inFC, ["SHAPE@"]) 
while lineRadians < 2*math.pi:
    polarCoords = PolarToCartesian(lineDistance, lineRadians)
    print "x = ", polarCoords[0] + orginX
    print "y = ", polarCoords[1] + orginY
    array = arcpy.Array([arcpy.Point(orginX, orginY),arcpy.Point(polarCoords[0] + orginX, polarCoords[1] + orginY)])
    polyline = arcpy.Polyline(array)
    cursor.insertRow([polyline])
    lineRadians = lineRadians + radiansBetweenLines
del cursor