r.klingeresri-de-esridist

The Heartbeat of a Region

Blog Post created by r.klingeresri-de-esridist Employee on Dec 28, 2017

Let's start with the result:

The heartbeat of BerlinOne of the major questions of a location is: What/Whome can I reach in [x] minutes? Therefore ArcGIS offers the Network Analyst with its tool "Service Areas". We can analyse the area reachable by car (other modes are possible) in a given time period. By using several input points we can also analyze patterns in a whole area.

Using a normal network would not alter the service areas during the day but by using the ArcGIS Online network dataset based on HERE data we can also use the traffi patterns for a day of the week and a given time. Traffi patterns can be analyswed in a 15min period. Analyzing the service areas during a day provides a valuable insight into the traffic situation of a region: big differences indicate first of all higher traffic volume on streets. The total area also indicates the quality of the network: large areas indicate a better accessibility using the car (maybe caused by a dense network of high-speed streets).

Automated Analysis

As I love geoprocessing, I was interested in the automation of the analysis and the visualisation as an animation. The process is therfore a two-stepped one.

First I created a simple point feature class called "StartPoints" which holds the locations of the "facilities" and added them to the Service Area Analysis Layer:

#create Analysis Layer
arcpy.na.MakeServiceAreaAnalysisLayer("https://www.arcgis.com/", "BeatCity Service Layer", "Driving Time", "FROM_FACILITIES", "5;10;15;20;25;30", "05.01.1900 00:01:00", "LOCAL_TIME_AT_LOCATIONS", "POLYGONS", "STANDARD", "DISSOLVE", "RINGS", "100 Meters", None, None)
#add locations from a point feature class
arcpy.na.AddLocations("BeatCity Service Layer", "Facilities", "StartPoints", "Name Description #;CurbApproach # 0;Attr_Minutes # 0;Attr_TravelTime # 0;Attr_Miles # 0;Attr_Kilometers # 0;Attr_TimeAt1KPH # 0;Attr_WalkTime # 0;Attr_TruckMinutes # 0;Attr_TruckTravelTime # 0;Breaks_Minutes # #;Breaks_TravelTime # #;Breaks_Miles # #;Breaks_Kilometers # #;Breaks_TimeAt1KPH # #;Breaks_WalkTime # #;Breaks_TruckMinutes # #;Breaks_TruckTravelTime # #", "5000 Meters", None, None, "MATCH_TO_CLOSEST", "APPEND", "SNAP", "5 Meters", "EXCLUDE", None)

Now we need to alter the times automatically to refelct the times of day. The start time is '05.01.1900 00:00:00' which is a Friday and the time is midnight.

We can access the properties quite easily. As I am working with a 3D scene this is my only map object in my ArcGIS Pro project:

#get the project
doc = arcpy.mp.ArcGISProject('current')
#get the map
map = doc.listMaps()[0]
#get the service area layer
sa_layer = map.listLayers("BeatCity Service Layer")[0]
#Get the solver properties object from the service area layer
solver_props = arcpy.na.GetSolverProperties(sa_layer)

As we now have access to the properties we can easily iterate and solve the service layer. The result will be exported as a 3D feature class:

import os
for hour in range(0,24):
     #create the datetime object
     date = datetime.datetime(1900, 01, 05 , hour,0,0)
     #we will use the string-version later on
     datestring = date.strftime("%d.%m.%Y %H:%M")
     #set the time of the day and the date:
     solver_props.timeOfDay=date
     #solve the network layer
     arcpy.na.Solve("BeatCity Service Layer", "SKIP", "TERMINATE", None, None)
     arcpy.AddMessage("copying hour" + str(hour))
     #as we would like to visualize this in 3D
     arcpy.ddd.FeatureTo3DByAttribute(r"BeatCity Service Layer\Polygons", arcpy.env.workspace + os.sep + "Friday_" + str(hour), "ToBreak", "ToBreak")
     #and add a field with the date for later usage if wanted:
     arcpy.AddField_management(arcpy.env.workspace + os.sep + "Friday_" + str(hour), "datetime","Date","", "",8, "datetime", "NULLABLE", "REQUIRED")
     arcpy.management.CalculateField(arcpy.env.workspace + os.sep + "Friday_" + str(hour), "datetime", '"' + datestring + '"', "PYTHON_9.3")

Warning: Solving the layer consumes credits!

The results are stored in seperate feature classes which is not optimal yet this is my solution at the very moment. For

Visualization and Export

As we have seperate layers now I am applying the same style to all of them and triggering the visibility of each layer and exporting the layout which needs to be created prior the export!

My Layout looks something like this:

layout with tilted view and extrusioned layer

I added a text symbol with a default Text "Friday" in the layout to export the timestamp in the png result as well. Furthermore I added a marker symbol with the cities name to the layout and tilted the view a bit to increase the depth effect for the extruded layer. But first we need to disable all layers for our automated export with a quite unsophisticated approach:

#deselcting all layers:
for layer in map_obj.listLayers():
     if layer .name[0] == "F":
          layer .visible = False

Now we can iterate through the layers, apply the unique value renderer, alter the symbols to get polygons without borders, set the extrusion, alter the text element to get the right timestamp and export the layout as png:

for hour in range(0,24):
     poly_layer = map_obj.listLayers("Friday_" + str(hour))[0]
     poly_layer.visible = True
     symbology = poly_layer.symbology
     symbology.updateRenderer('UniqueValueRenderer')
     symbology.renderer.fields = ['ToBreak']
     symbology.renderer.colorRamp = doc.listColorRamps('Green-Blue (6 Classes)')[0]
     #alter all symbols:
     for sym in symbology.renderer.groups[0].items:
          sym.symbol.outlineColor = {'RGB': [255, 255, 255, 0]}
          sym.symbol.size = 0.0
     poly_layer.symbology = symbology
     #extrude:
     poly_layer.extrusion('BASE_HEIGHT', "3000-100*[ToBreak]")
     #export:
     datestring = "2017/12/08 " + str(hour) + ":00"
     for lyt in doc.listLayouts():
          for elm in lyt.listElements("TEXT_ELEMENT"):
               if elm.text[0] == "2":
                    elm.text = datestring
     lyt.exportToPNG(arcpy.env.workspace + os.sep + "friday_seattle_final" + str(hour) + ".png", 300)
     poly_layer.visible = False

Unfortunately the whole rendering and exporting process takes about 2.5minutes per scene on my Lenovo X270 with no dedicated GPU.

The result is a bunch of pngs that can be converted intzo a gif using several gif-creation  sites on the web or the open source image editor GIMP.

The result looks like this:

 

You can also append all polygon features from the result, add the datetime attribute, fill the used time, time-enable the feature class, and use the animation options in ArcGIS Pro to export the whole animation as a GIF:

The result looks quite similar to the file-wise export. Yet you don't have all the possibilities of the layout view. But still the result is quite good (added 5, 10, 15, 20 mil buffers):

Outcomes