ArcGIS Pro: A year later

I mentioned previously that I was finally considering making the jump from ArcMap to ArcGIS Pro. I made the switch in early 2018, and after an adjustment period of a few weeks, I have not looked back.

Seriously.

For those still waiting to make the switch, just do it. I highly recommend completing the official tutorials. It is a huge adjustment to go from dialogue boxes to the ribbon interface, akin to the user shock experience when switching over from Office 2003 to Office 2007.

You’ll also need to learn how to deal with multiple data  map frames and layouts in one map document project, the lack of a standalone ArcCatalog, and the ability to link 2D and 3D views in one place. There are also new terms which you will need to learn, and new concepts.

A few of my favourite things:

• Seamless integration with AGOL/Portal
• Labelling properties now accessible through groups on the ribbon – no more Inception levels of dialogue boxes
• Data driven pages is more robust
• BIM integration
• arcpy.mapping overhauled to arcpy.mp – a much more intuitive module
• One ArcGIS project per client containing multiple data frames for each task – no need for dozens of map documents
• Python 3 support (and built in Anaconda!)

I could go on. Our local licence server for Desktop went down a few weeks ago, and I only found out a few days after it happened because I just don’t open ArcMap anymore. With the latest update to Pro, I can now quickly launch a project without needing to save it first, a feature that was lacking up until now. My last reason to use ArcMap is gone.

A rant about attachments

I remember when attachments were first introduced in ArcGIS Desktop (10.2 I think? whoops it was ArcGIS 10). It was a very useful feature, and more functionality was added over the years.

It also made mobile data capture even easier. The fieldworkers would go out, do their assessments, and attach multiple photos to their points. However, attachments with Collector has caused me so much frustration. Specifically, syncing with attachments.

The nature of the work we do (and the economic environment we are in) means that by default, I take the maps offline so that the fieldworkers can carry out their assessments, and then sync back to AGOL when they are on lunch break (or whenever they can pick up WiFi). I discovered a few years ago that once one hits a certain threshold (like 20 attachments in the map), there are going to be problems syncing.

It will just outright fail, or take very long and may need to be attempted a number of times. Why is this? I don’t know. Over the years, I’ve encountered this issue on all types of devices – the latest iPhones, low-end Android tablets, high-end Android tablets, mid-range Android phones…

What it seems like to me is that Collector “expects” a certain connection speed, and when it doesn’t get it, it times out and rolls back the sync. Fair enough – I’ve found multiple delta tables on devices I’ve needed to recover the databases from due to failed sync attempts. On a current project, they are using rugged devices which have really awful network chips (as in, I need to stand about 1 or 2m away from the access point so that I can take the maps offline). Naturally, at the end of the first day, each device had dozens of features with multiple attachments each, which refused to sync.

They have been out in the field for 2 weeks. Everyday, I have to manually retrieve the databases from the device, recover them, and push them out into appropriate geodatabases once I’ve determined what’s inside them.

So clear

I can deal with all of that, because Python is a tool that I maaay have mentioned here before. What I cannot deal with is the fact that attachments are still lost during geoprocessing. The fact that it was added as an environment setting in ArcGIS 10.5 and has been available in ArcGIS Pro for a while is of little comfort to me as I currently have access to neither.

Fine. I store the GlobalIDs in another field, merge the features together into their correct feature classes, enable attachments and insert the records from the corresponding attachment tables. Of course, I forget that the relationship class is now messed up, as it’s linking through the (now incorrect) GlobalID fields instead of the fields I stored the original IDs in.

After staring at the screen cross-eyed, I then realise that I only need to provide the attachments as jpgs in a folder, which I can extract from the tables using the original IDs and write into subfolders based on the feature type. I don’t actually need to link them back together since the technician does not need to view the photos to complete the work in ArcMap. /endrant

Create points without XY from a table in ArcPy

I had a point feature class containing features with unique names. I also had a table containing hundreds of records, each “linking” to the point feature classes via name. Normally, I could create points for the records by using Make Query Table.

I say “linking” because from looking at the data, I could see which records belonged to which point, but unfortunately there were spelling mistakes and different naming conventions e.g. if the point was called “ABC Sewage Treatment Works”, some of its matching records in the table would be “ABC WWTP”, “AB-C Waste Water Treatment Plant”.

	'''
	@author Cindy Jayakumar
	@date 31/01/2017
	– Inserts a new point with the selected geometry in to_lyr
	– Adds attributes to that point from from_lyr
	– Updates the record in from_lyr
	Uses selected features
	For the python window in ArcMap
	'''
	import arcpy
	from_lyr = r'C:\Some\Arb\Folder\work.gdb\tbl_test'
	mxd = arcpy.mapping.MapDocument("CURRENT")
	out_lyr = arcpy.mapping.ListLayers(mxd, "lyr")[0]
	'''
	– Select the single point in to_lyr that will have its geometry duplicated
	– Select the rows in from_lyr that you want to be inserted into to_lyr
	'''
	def duplicateAssets(to_lyr):
	# Matching point is selected in to_lyr
	geom = [geo[0] for geo in arcpy.da.SearchCursor(to_lyr, "SHAPE@")][0]
	# Create insert cursor on same layer
	inscursor = arcpy.da.InsertCursor(to_lyr, ("SHAPE@", "Unique_ID"))
	with arcpy.da.UpdateCursor(from_lyr, ("Unique_ID", "Matched")) as cursor:
	for row in cursor:
	# Build the new point to be inserted
	point = [geom, row[0]]
	inscursor.insertRow(point)
	print("Inserted " +str(int(row[0])))
	# Update the record with the name of the layer the point was inserted into
	row[1] = to_lyr.datasetName
	cursor.updateRow(row)
	del inscursor
	duplicatePoint(out_lyr)

view raw

duplicatePoint.py

hosted with ❤ by GitHub

Add a new field to a feature class using NumPy

I needed to add a field to dozens of layers. Some of the layers contained the field already, some of them contained a similar field, and some of them did not have the field. I did not want to batch Add Field, because not only would it fail on the layers which already had the field, but it is super slow and I would then still have to transfer the existing values from the old field to the new field.

	'''
	@date 12/10/2015
	@author Cindy Williams
	Adds a new field to layers in a map document, based
	on a current field.
	For use in the Python window in ArcMap.
	'''
	import arcpy
	import numpy
	mxd = arcpy.mapping.MapDocument("CURRENT")
	lyrs = arcpy.mapping.ListLayers(mxd)
	# Create a numpy array with a link oid field and the fields to be added
	narray = numpy.array([], numpy.dtype([('objid', numpy.int), ('GIS_ID', '|S10'),]))
	for lyr in lyrs:
	if lyr.isFeatureLayer:
	# Find the current field with the values
	field = [field.name for field in arcpy.ListFields(lyr, "GIS_*")][0]
	# Prevent the tool from failing
	if field != "GIS_ID":
	# Add the field
	arcpy.da.ExtendTable(lyr,"OID@", narray, "objid")
	# Copy the field values to the new field
	with arcpy.da.UpdateCursor(lyr, ("GIS_ID", field)) as cursor:
	for row in cursor:
	row[0] = row[1]
	cursor.updateRow(row)

view raw

addFieldWithExtendTable.py

hosted with ❤ by GitHub

I wrote the tool in the ArcMap Python window, because I found it easier to load my layers into an mxd first as they were lying all over the place. The new field to be added is set up as a NumPy array, with the relevant dtype.

The script loops over all the layers in the document, adding the field via the Extend Table tool, and then transferring the values from the old field to the new field. Deleting the old field at the end would be an appropriate step to include, but I didn’t, purely because I’ve lost data that way before.

Split line interactively using ArcPy FeatureSets

Back when I was doing a lot of GIS work for civil and electrical asset management, a common task I had was to ensure that all civil structures represented as a polyline (reticulation and bulk pipelines, eletrical feeder cables, road surfaces etc) were split at the physical road intersections as seen on the latest available aerial imagery, as well as where they crossed each other.

This task wasn’t too much trouble if the lines were already drawn across each other, as more often than not, the lines crossed at the road intersections. The result was usually achieved by using the Planarize option on the Advanced Editing toolbar. However, it became a problem for new areas which had to be digitised by hand (no CAD drawings or shapefiles received), or lines which had to be corrected.

I would have to visually check on the imagery where the road intersections were, and then manually split each polyline. I ran through some ideas, including digitising all road intersections as points (which would be good data to have anyway), and then splitting all the lines at the closest intersection by using the Split Line at Point tool.

I did however want to create a method of splitting the line multiple times on the fly, instead of just splitting it once (which is the option provided by the Split button on the Editor toolbar). I wanted to be able to select a line, click as many points along the line as I wanted, and then split the line into those points.

	'''
	@date 06/10/2015
	@author Cindy Williams
	@version 1.0
	Uses an arcpy FeatureSet to allow the user to interactively create
	points on a line which will be used to split that line into parts,
	while retaining all attributes.
	Requirements:
	– Works as a script tool. Add to a toolbox.
	– Can only work on a single line for now. Select the line, right
	click the layer and "Create layer from selected features". Use
	this selection layer as input to the tool.
	– Digitise points snapped to the selected line. Tool will fail if not
	snapped
	– Digitise points left to right along the line.
	Future versions will include:
	– split multiple lines at same time
	– snap points automatically/error that points should be snapped
	– points will be sorted after digitising to ensure correct order
	Built on an idea here http://gis.stackexchange.com/questions/154708/standalone-python-script-to-split-a-polyline-with-a-point-layer
	'''
	import arcpy
	inset = arcpy.GetParameter(0) # FeatureSet of interactive points from ArcMap
	inline = arcpy.GetParameterAsText(1) # Input polyline feature layer
	# Get all field names excluding dynamic fields
	fields = ["SHAPE@"]
	fields.extend([field.name for field in arcpy.ListFields(inline) if field.name not in ["OBJECTID", "Shape_Length", "Shape"]])
	# Open an insert cursor for inserting the cut lines
	cursor_ins = arcpy.da.InsertCursor(inline, fields)
	# Get the spatial reference from the polyline
	sr = arcpy.Describe(inline).spatialReference
	# Get the polyline and delete the original record
	with arcpy.da.UpdateCursor(inline, fields) as cursor:
	for row in cursor:
	polyline = row
	cursor.deleteRow()
	# Iterate over the points
	with arcpy.da.SearchCursor(inset, "SHAPE@") as cursor:
	for row in cursor:
	point = row[0].firstPoint
	# Build a temporary diagonal line from the point
	diag = arcpy.Polyline(arcpy.Array([arcpy.Point(point.X + 10.0, point.Y + 10.0), arcpy.Point(point.X-10.0, point.Y-10.0)]), sr)
	# Cut the original polyline
	geom = polyline[0].cut(diag)
	# Create a new line with the right side of the geom i.e. finished with this side
	newline = [geom[1]]
	newline.extend(polyline[1:])
	# Insert the line into the feature layer
	cursor_ins.insertRow(newline)
	# Set the left side of the cut line as the new line for cutting
	polyline = [geom[0]]
	polyline.extend(newline[1:])
	# Insert the last portion of the line
	cursor_ins.insertRow(polyline)
	del cursor_ins
	del polyline

view raw

splitLineAtInteractivePoint.py

hosted with ❤ by GitHub

The tool iterates over all the points created by the user in the FeatureSet. For each point, a temporary diagonal is created to cut the original polyline. This new part is saved in a list and the tool runs again on the original polyline (minus the new part). Once all the points have been processed, the original polyline is removed from the feature class, and the new parts are inserted. All feature attributes are maintained.

My intention with this tool was to expand it with more features, fix the bugs (it would randomly not work) and upload it as a toolbox. However, the need to do that fell away. I also realised I was starting to get into a space where I was forcing as much code as I could into my work in order to get through the day.

Convert spreadsheets to GDB tables and display XY

I’m sick today, so I decided to choose one of the simpler scripts for today’s post. This one automates a common task of the GIS professional: converting a spreadsheet of coordinates to GIS. The script assumes you’ve formatted each sheet appropriately.

	#
	# @date 28/07/2015
	# @author Cindy Williams
	#
	# Converts all the sheets in a spreadsheet
	# to file gdb tables and displays the XY
	# coordinates in ArcMap.
	#
	# Assumes the spreadsheet has been formatted
	# according to GIS standards.
	#
	# For use in the Python window in ArcMap.
	#
	import arcpy
	import xlrd
	arcpy.env.workspace = r"C:\Some\Arb\Folder\work.gdb"
	sr = arcpy.SpatialReference(4326)
	xls = r"C:\Some\Arb\Folder\Coords.xlsx"
	wb = xlrd.open_workbook(xls)
	sheets= [sheet.name for sheet in wb.sheets()]
	field_x = "Long"
	field_y = "Lat"
	for sheet in sheets:
	print("Processing " + sheet)
	tbl = "tbl_" + sheet
	try:
	arcpy.conversion.ExcelToTable(xls, tbl, sheet)
	print("Converted " + sheet)
	arcpy.management.MakeXYEventLayer(tbl, field_x, field_y, sheet, sr)
	print("Displaying XY")
	except Exception as e:
	print(e)
	print("Script complete.")

view raw

convertSheetsAndDisplayXY.py

hosted with ❤ by GitHub

For each sheet in the workbook, the script creates a file geodatabase table, and displays the coordinates as a layer in ArcMap using the WGS 1984 spatial reference. The script is very basic, and therefore will be easy to change to suit your own data.

Export selected data driven pages only

A common task I perform with the asset management data is to create mapbooks per service (Water, Electricity etc) per main town/area/grid. I normally do this at the start of the project, to display the data we currently have. The asset guys then take the maps to the client, who will mark up amendments and additions.

While I cringe at having to create hard copies of anything, it’s the quickest way to get the information out of the client. It’s also relatively easy for me to create the mapbooks – I just have to open my service templates, save a new copy, set the data sources and CS to my current area, and export. I could automate this part of the process also, but often I need to visually confirm that the “current” data I am using is actually “current”.

When I set up the data driven pages for an area, I switch on all the asset data layers I have, and see which towns they are grouped into. I create best-fit rectangles around them, either by eyeballing it and drawing it in graphics, or using the Strip Map Index Features tool. Often, some areas won’t have Electricity data, or the Water Supply data only covers 3 out of the 5 grid blocks which covers the town.

Before, I would manually check this, and create separate index layers for each service, which would be used for that service’s data driven pages. This was a terrible way to do things, but doing it that way is what made me determined to find a better way. To Python!

	#
	# @date 13/07/2015
	# @author Cindy Williams
	#
	# Intersects feature layers with the index layer
	# for selected data driven export to a multi-page
	# PDF with embedded layers.
	#
	# For use as a standalone script.
	#
	import arcpy
	import os
	map_folder = r"C:\Some\Arb\Folder\mxd
	pdf_folder = r"C:\Some\Arb\Folder\pdf"
	kml_folder = r"C:\Some\Arb\Folder\kml"
	# Get all the mxds in the map folder
	(_,_, mapdocs) = os.walk(map_folder).next()
	for md in mapdocs:
	mxd = arcpy.mapping.MapDocument(os.path.join(map_folder, md))
	df = arcpy.mapping.ListDataFrames(mxd)[0] # First data frame
	ddp = mxd.dataDrivenPages
	ddp_lyr = ddp.indexLayer
	# Only get the layers in the Data group
	data_lyrs = [lyr for lyr in arcpy.mapping.ListLayers(mxd) if lyr.isFeatureLayer and lyr.longName.split("\\")[0] == "Data"]
	for lyr in data_lyrs:
	# Selects the data driven pages that intersects all feature layers
	# in the Data group
	arcpy.management.SelectLayerByLocation(in_layer=ddp_lyr,
	overlap_type="CONTAINS",
	select_features=lyr,
	selection_type="ADD_TO_SELECTION")
	print("Data driven pages selected: {}.".format(len(ddp.selectedPages)))
	pdf_name = os.path.join(pdf_folder, md[:-4] + ".pdf")
	# Uses the PDF export function from the DataDrivenPagesClass
	ddp.exportToPDF(out_pdf=pdf_name,
	page_range_type="SELECTED",
	multiple_files="PDF_SINGLE_FILE",
	resolution=100,
	layers_attributes="LAYERS_AND_ATTRIBUTES",
	georef_info="True")
	print("Exported " + md)
	# Cannot export to KML with basemap in mxd. Remove layer
	# first and add back in
	kmz_name = os.path.join(kml_folder, md[:-4] + ".kmz")
	# Converts the mxd to KML. Only the Data group layers are
	# switched on, so only they will be exported. Can enforce
	# the layer visibility by lyr.visible = True
	arcpy.conversion.MapToKML(in_map_document=mxd,
	data_frame=df.name,
	out_kmz_file=kmz_name,
	map_output_scale=10000)
	print("KML conversion complete.")
	print("Script complete.")

view raw

selectVisibleFeaturesForDDPExport.py

hosted with ❤ by GitHub

This script takes a folder containing the various service MXDs, which have been set up so that the operational data is in a group layer called Data. Each of these layers is intersected with the data driven pages index layer, to determine which DDP polys contain the operational data. In Line 35, I’ve specified the selection_type as ADD_TO_SELECTION, so as not to lose the previously selected pages.

In Line 39, instead of making the traditional call to arcpy.mapping.ExportToPDF, I used dataDrivenPages.exportToPDF instead. The funny thing about the ArcGIS documentation is that I can look at the same pages so many times, and still discover new things.

What I like about this method is that it allows a SELECTED page range to be exported. This allows me to only have one index layer for an area 🙂 At the end of the script, I put some code in to convert the entire map to a KMZ file. That portion currently does not work, because the imagery basemap needs to be removed from the mxd before it can be converted. This should be as simple as calling RemoveLayer on the basemap, but I always have fantastic problems working with the layer movement methods (insert, add, move etc).

Edit: I just saw that I actually posted this script in July as well, when I wrote it initially. I was just so excited.

Using a query table to represent a 1:M relationship spatially

I first discovered (and used) the Make Query Table tool during my second week at Aurecon (March 13 2012, according to my OneNote). This was about a month before I started using ModelBuilder (to combat the frustrations of ArcMap), and about 6 months before I started Python (to combat the frustrations of ModelBuilder).

I’m just giving a bit of context because this was before everything clicked into place for me. Before this point, I treated everything I learnt in my studies and during my internship as separate silos of information: GIS, Databases, Programming.

I did not even realise until after I was working at Aurecon that the query expressions used in ArcMap are SQL, despite me studying all those things. It just shows how one’s mindset can block progress, and how I allowed the awful experience of learning Computer Science in Afrikaans to stop me from letting things “click” for so long.

Back to the tool. Joins in ArcGIS are notoriously slow, and 1:M joins are not allowed (technically they are, but in the sense that an arbitrary matching feature will be joined). Naturally, the relationship between the GIS and the asset register is 1:M.

For example, a single point is used to represent the physical, spatial location of an asset, reservoir WRV-00001. In the asset register, this reservoir is unbundled into various components – storage tank, building, fence etc. Each of these assets have their own unique asset ID, but all have the same GIS ID.

I now need to represent all these assets spatially. The points/lines will all be on top of each other, but that’s fine. The Make Query Table tool does exactly this, but it is…quirky. I’ve compiled a list of things to remember when using this tool (supplemented by this question on GIS.SE):

1. Tables/feature classes in the relationship should be stored within the same database: I tend to remember this step only after I add the inputs and the tool shouts at me.
2. Add the feature class first in the multivalue parameter: The format of the input relies on the format of the first argument in the multivalue parameter control. The feature class should be added first to ensure that the output is a layer, otherwise it will be a table view.
3. Enclose table names and field names in quotes: For example, I wish to join the asset register table ar to the point layer points using their common field GIS_ID. By default, the tool encloses my whole expression in quotes "points.GIS_ID" = "ar.GIS_ID". This will cause the tool to fail. Add extra quotes around the field names "points"."GIS_ID" = "ar"."GIS_ID".
4. Choose the NO_KEY_FIELD option: Trying to add key fields causes some erratic behaviour (???). Just don’t do it. By selecting this option, the existing ObjectID field from the input will be used.
5. The output layer will appear to have no symbology: Go into the Layer properties, click the Symbology tab, click the existing symbol, then OK, OK. It’s a bug./li>
6. Persist to disk!: Remember to export the layer to a feature class, otherwise the layer will only exist in the map document.

ArcMap Woes Part 6: Random failure of gp tools on grouped layers

I had a series of errors recently while I was trying to run a few geoprocessing tools manually on layers within a group layer. The error message kept saying it couldn’t open the layer, even though the layer was clearly in the map. Even restarting ArcMap did not help.

I eventually discovered that geoprocessing tools may randomly fail on layers which are part of a group layer. I say randomly because sometimes the tools would run as expected, other times not. I’ve only experienced this anomaly when using the ArcMap interface – running the same sequence of tools on all the layers in a group layer in a for loop in Python does not give this error.

ArcMap Woes Part 5: Exporting map without options

Sometimes, when in ArcMap, and one wants to export a map as a JPG or PNG or something, there won’t be any options. You can’t change the resolution, can’t change the export quality, can’t include a world file or anything. This used to frustrate me enormously, because then I would have to restart ArcMap to get the options to appear.

Sometimes even that wouldn’t work. I would then have to disturb a colleague and have them open the same mxd, to find that the options have miraculously reappeared. I have since figured out that if this happens, one needs to cancel the export dialog, switch to layout view and then export the map.

I am aware that it is possible to export the map with options while in data view; it’s just a random thing that happens. At least I know how to solve this particular issue.