Sort a polygon feature class spatially using Python

I was manually adjusting a custom grid for a map series this morning, a task which I have not had to do for several years now, when the time came to export the maps. By this point, my grid numbers were completely messed up because of all the shifting and reshaping I had to do.

It would make sense to number the grid from left to right, top to bottom. However, the ability to spatially sort a feature class in ArcGIS is hidden behind an Advanced Licence in the Sort tool.

License: For the Field(s) parameter, sorting by the Shape field or by multiple fields is only available with an Desktop Advanced license. Sorting by any single attribute field (excluding Shape) is available at all license levels.

I had this prickly feeling in the back of my head that I had found a workaround for this before. I went searching through the blog (since I have posted about things twice because I forgot about the first post), but I didn’t find it.

Undeterred, I searched through my gists (unrelated note – why does using the “user:<username> <search term> ” trick in the GitHub search bar now flag as spam?!) and I found it!

	#
	# @date 13/07/2015
	# @author Cindy Williams
	#
	# Sorts a polygon feature class by descending X and ascending Y,
	# to obtain the spatial order of the polys from left to right,
	# top to bottom. Bypasses the Advanced licence needed to
	# do this in the Sort tool.
	#
	# For use in the Python window in ArcMap.
	#
	import arcpy
	mxd = arcpy.mapping.MapDocument("CURRENT")
	lyr = arcpy.mapping.ListLayers(mxd, "DDP")[0]
	# Get the extent of each feature, and store as tuple (name, x coord of top right, y coord of top right)
	ddp = [(row[0], row[1].extent.XMax, row[1].extent.YMax) for row in arcpy.da.SearchCursor(lyr, ("Name", "SHAPE@"))]
	# Sort by x descending
	ddp.sort(key=lambda row:row[1], reverse=True)
	# Sort by y ascending
	ddp.sort(key=lambda row:row[2])
	# Reverse the list to get the L->R, T->B order
	dct = {}
	# Must still be optimised
	for i, j in enumerate(reversed(ddp)):
	# Key is name, value is index + 1
	dct[j[0]] = i + 1
	# Write values back to feature class
	with arcpy.da.UpdateCursor(lyr, ("Name", "MapNum")) as cursor:
	for row in cursor:
	row[1] = dct[row[0]]
	cursor.updateRow(row)

view raw

sortFeatureClassByExtent.py

hosted with ❤ by GitHub

Link for mobile users is here.

I’m pretty sure there is a much better way to do this, but I need this so rarely that I don’t see the need to rewrite it. It also would work much better as a function. It takes a polygon feature class and sorts the selected features from left to right, top to bottom.

Description of this hastily drawn, blurry PowerPoint graphic: Chaotic, unsorted grid on the left. Slightly less chaotic, sorted grid on the right.

Inserting spaces into CamelCase

I have some trauma related to regular expressions from my university days. Yes, they are super useful. Do I enjoy using them? Does anyone?

I routinely modify and create new instances of the AMIS GIS data model. Depending on the client and the type of assets they have, there can be 4 feature datasets, containing about 5 feature classes each, or 9 feature datasets with up to 50 feature classes spread across the database.

Inevitably at some point in the process of losing my mind, I will forget that the alias is lost when creating a new feature class, or randomly copying it over or whatever. I then end up with dozens of feature classes with really ugly looking abbreviated layer names in CamelCase. So for posterity, and so that I will never ever forget this ridiculously easy thing again:

	#
	# @date 05/10/2017
	# @author Cindy Jayakumar
	#
	# Change the alias of feature classes to the
	# title case version of the CamelCase name
	#
	# e.g. WaterPumpStation => Water Pump Station
	import arcpy
	import re
	gdb = r"C:\Some\Arb\Folder\test.gdb"
	def alterFCAlias(name):
	return re.sub(r"\B([A-Z])", r" \1", name)
	for root, _, fcs in arcpy.da.Walk(gdb):
	for fc in fcs:
	alias = alterFCAlias(fc)
	arcpy.AlterAliasName(fc, alias )
	print("{0} => {1}".format(fc, alias))

view raw

alterCamelCaseAlias.py

hosted with ❤ by GitHub

	# @credit https://stackoverflow.com/a/199215
	import re
	text = "ThisIsATest"
	re.sub(r"\B([A-Z])", r" \1", text)

view raw

insertSpaceBetweenCamelCase.py

hosted with ❤ by GitHub

One thing I’ve omitted from the script is that I usually store the feature class names with a prefix, e.g. fc_name = wps_WaterPumpStation. In this case, I would use split on the feature class name before passing it to the alterFCAlias function i.e. fc_name.split("_")[1].

Link for mobile users is here.

Aaaaaaand I’ve just realised that I’ve forgotten this basic task so many times over the last few years that I actually already have a blog post about it, except for some reason I was updating the layer names in ArcMap every time instead of resetting it once on the feature classes themselves.

Convert a list of field names and aliases from Excel to table using ArcPy

I went digging through my old workspace and started looking at some of my old scripts. My style of coding back then is almost embarrassing now 🙂 but that’s just the process of learning. I decided to post this script I wrote just before ArcGIS released their Excel toolset in 10.2.

	'''
	@date 09/07/2013
	@author Cindy Williams
	Converts a spreadsheet containing field names and aliases
	into a file geodatabase table.
	'''
	import arcpy
	xls = r"C:\Some\Arb\Folder\test.xlsx\Fields$"
	tbl = r"C:\Some\Arb\Folder\test.gdb\tbl_Samples"
	with arcpy.da.SearchCursor(xls,("FIELD_NAME","ALIAS")) as cursor:
	for row in cursor:
	arcpy.management.AddField(tbl, row[0], "DOUBLE", field_alias=row[1])
	print("Adding field {}…".format(row[0]))
	print("Fields added successfully.")

view raw

createTableFieldsFromXLS.py

hosted with ❤ by GitHub

From what I can recall, I needed to create a file geodatabase table to store records of microbial sample data. Many of the field names were the chemical compound themselves, such as phosporus or nitrogen, or bacterial names. For brevity’s sake, I had to use the shortest field names possible while still retaining the full meaning.

I set up a spreadsheet containing the full list of field names in column FIELD_NAMES and their aliases in ALIAS. I created an empty table in a file gdb, and used a SearchCursor on the spreadsheet to create the fields and fill in their aliases.

This solution worked for me at the time, but of course there are now better ways to do this.

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.

Write feature class properties to a Word document using python-docx

It’s been too long since I’ve posted some code. A few months ago I had a requirement to basically create a UML view of a file gdb and present each feature class as tables for inclusion in a functional specification Word document.

Sadly, since the demise of ArcGIS Diagrammer, there has never been something to take its place. The actual request I received was to take screenshots of the properties dialog box of each feature class in ArcCatalog, and paste those into the document with appropriate headings.

I recognised this request for the total waste of my time it would be, and promptly set about looking for an alternative. I realised I had python-docx already installed.

	'''
	@date 15/09/2015
	@author Cindy Williams
	Uses the python-docx package to write various
	properties of a feature class to docx format.
	'''
	import arcpy
	from docx import Document
	arcpy.env.workspace = r"C:\Some\Arb\Folder\work.gdb"
	doc = r"C:\Some\Arb\Folder\test.docx"
	document = Document()
	document.add_paragraph("This document contains a description of fields per feature class.")
	for fc in arcpy.ListFeatureClasses():
	p = document.add_paragraph("Feature class: ")
	p.add_run(fc).bold = True
	p.add_run("\n")
	table = document.add_table(rows=1, cols=2, style='Table Grid')
	header_cells = table.rows[0].cells
	header_cells[0].text = "Field Name"
	header_cells[1].text = "Field Type"
	for field in arcpy.ListFields(fc):
	row_cells = table.add_row().cells
	row_cells[0].text = field.name
	row_cells[1].text = field.type
	document.add_page_break() # Each feature class on its own page
	document.save(doc)

view raw

writeFieldsToDOC.py

hosted with ❤ by GitHub

For each feature class in the geodatabase, a new paragraph is started with the name of the feature class in bold. A line break is inserted, followed by a table. The name and type of each field is added as a new row into the table, and a page break is inserted to start the properties of the next feature class on a new page.

It took me about an hour to look up alternative methods and to put this script together. Most of the time was spent on getting the cells in the table to insert properly, using the age-old method of trial and error. I did it this way to save myself the pain of manually inserting screenshots, knowing that if the format of the feature classes changed I would have to do new screenshots repeatedly.

Display raster properties arcpy vs numpy

	#
	# @date 25/08/2015
	# @author Cindy Williams
	#
	# Prints a comparison of the specified raster
	# properties using arcpy and numpy.
	#
	# For use as a standalone script.
	#
	import arcpy as ap
	import numpy as np
	ext_lyr = ap.management.MakeFeatureLayer(r"") # Extent polygon shapefile/feature class
	input_raster = ap.Raster(r"") # Raster containing values
	def get_raster_properties(ras, prop):
	array = ap.RasterToNumPyArray(ras, nodata_to_value=0)
	ap_answer = ap.management.GetRasterProperties(ras, prop)
	if prop == "MEAN":
	np_answer = array.mean()
	if prop == "SUM":
	np_answer = array.sum()
	if prop == "ROWCOUNT":
	np_answer = array.shape[0]
	if prop == "COLUMNCOUNT":
	np_answer = array.shape[1]
	# etc
	return (ap_answer, np_answer)
	with arcpy.da.SearchCursor(ext_lyr, ("SHAPE@", "Name")) as cursor:
	for row in cursor:
	# Get the extent of the current feature
	ext = row[0].extent
	# Build the feature envelope from the extent
	enve = "{0} {1} {2} {3}".format(ext.XMin, ext.YMin, ext.XMax, ext.YMax)
	# Clip the raster to the envelope in memory
	clipped_raster = ap.Raster(ap.management.Clip(input_raster, enve, r"in_memory\clip"))
	# Get the properties
	(arcpy_ans, numpy_ans) = get_raster_properties(clipped_raster, "MEAN")
	print("{0}\tarcpy: {1}\tnumpy: {2}".format(row[1], arcpy_ans, numpy_ans))
	# Clean up after each run
	ap.management.Delete("in_memory")

view raw

displayRasterProperties.py

hosted with ❤ by GitHub

Pan and zoom all data frames to match the data driven data frame in ArcMap

A colleague had a map document containing 9 data frames showing various climate scenarios. He asked me for a script to pan and zoom each data frame to the extent and scale of the data driven data frame. The data driven pages were set up beforehand and must be enabled manually.

	#
	# @date 25/08/2015
	# @author Cindy Williams
	#
	# Pans and zooms each data frame in an ArcMap document
	# to the extent and scale of the data frame containing
	# the data driven pages index layer before exporting.
	#
	# For use as a standalone script.
	#
	import arcpy
	import os
	map_doc = r"" # Your mxd here
	mxd = arcpy.mapping.MapDocument(map_doc)
	out_folder = r"" # Your folder here
	if mxd.isDDPenabled: # Check if data driven pages are enabled in the map
	ddp = mxd.dataDrivenPages
	ddp_df = ddp.dataFrame # Get the data frame containing the data driven pages index layer
	# Get all data frames besides the one containing the index layer
	dfs = [df for df in arcpy.mapping.ListDataFrames(mxd) if not df.name == ddp_df.name]
	# Loop over the data driven pages
	for i in range(1, ddp.pageCount + 1):
	ddp.currentPageID = i
	cur_ddp = ddp.pageRow.getValue(ddp_name)
	ddp_jpeg = os.path.join(out_folder, str(cur_ddp) + ".jpg")
	for df in dfs:
	df.extent = ddp_df.extent
	df.scale = ddp_df.scale
	arcpy.mapping.ExportToJPEG(mxd, ddp_jpeg, resolution=200)
	print("Exported " + ddp_jpeg)
	else:
	print("Please enable Data Driven Pages before continuing.")

view raw

panZoomDataFramesToDDPDataFrame.py

hosted with ❤ by GitHub

In Line 21, the data frame containing the data driven pages index layer is stored, so that it can be excluded from the list of all the data frames in the document in Line 24. When looping over the data driven pages for export, the scale and extent of each data frame are set to match those of the main data frame before writing to jpeg.

Convert multi-page PDF file to TIFF

	#
	# @date 18/08/2015
	# @author Cindy Williams
	#
	# Converts a range of pages from a pdf to separate tiff
	# files using the PDF to TIFF tool added at ArcGIS 10.3.
	# For very basic conversion, with minimal error handling.
	#
	# For use as a standalone script.
	#
	import arcpy
	import os
	def convertPDFtoTIFF(in_pdf, out_folder, out_tiff, start_page=1, end_page=1):
	if start_page > end_page:
	print("The start page number cannot exceed the end page number.")
	else:
	for i in range(start_page, end_page + 1):
	out_tiff_name = os.path.join(out_folder, out_tiff + str(i) + ".tif)
	arcpy.conversion.PDFtoTIFF(in_pdf, out_tiff_name, pdf_page_number=i)
	print("Processed page " + str(i))
	pdf_in = r"C:\Some\Arb\Folder\test.pdf"
	tif_folder = r"C:\Some\Arb\Folder\tiffs"
	convertPDFtoTIFF(pdf_in,tif_folder, "Page_", 2,9)
	print("Script complete.")

view raw

convertMultiPagePDFtoTIFF.py

hosted with ❤ by GitHub

A new tool was added to ArcGIS 10.3 to convert a page in a PDF to TIFF. This is especially useful for the asset management work I do, as I often get pdf drawings or maps which need to be georeferenced. I created a function to wrap around the ArcGIS one to enable the conversion of more than one page in the pdf at a time.

By default, the function will only convert the first page. Specifying a start and/or end page will override the defaults. Link to the code is here.

Reformat all layer names in an ArcMap document

*Before I get into today’s script, I’ve noticed that the embedded gists do not appear on mobile devices when accessing the feed for this blog via Feedly or some other RSS reader. The post needs to be opened in a web browser to view the code. I’ll also be including a link to it at the end of the post.

When I create feature classes for asset management data manually instead of automatically, I always make sure to name them consistently. However, I rarely change the feature class alias manually. This means that when I need to pull all the layers into a map to create a kml, I need to rename each layer.

Which is why I turn to Python instead.

	#
	# @date 14/08/2015
	# @author Cindy Williams
	#
	# Quickly formats the names of all the layers in a map
	# document. In this example, the part of the name before
	# the first full stop is removed, underscores are replaced
	# spaces, a regex is used to insert spaces before capitals
	# and the string is converted to proper case.
	#
	# For use in the Python window in ArcMap.
	#
	import arcpy
	import re
	mxd = arcpy.mapping.MapDocument("CURRENT")
	lyrs = arcpy.mapping.ListLayers(mxd)
	for lyr in lyrs:
	new_name_replaced = lyr.name.split(".")[1].replace("_", " ")
	new_name_spaced = re.sub(r"(\w)([A-Z])", r"\1 \2", new_name_replaced)
	lyr.name = new_name_spaced.title()
	arcpy.RefreshTOC()

view raw

reformatAllLayerNames.py

hosted with ❤ by GitHub

Of course, the way the name gets split will change depending on the current format of the names. In this instance, the part of the name that I needed occurred after the first full stop, and contained underscores which needed to be replaced by spaces. I then used a regular expression to insert spaces between InitialCaps words.

Link to the code is here.