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.

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.

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

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 4: Permissions errors

Sometimes when I’m in ArcMap, I can’t start an edit session. Normally, I restart ArcMap and everything is fine.

Recently, for about a week, I could not edit inside any file gdb, or create new feature classes inside of an existing one or a new one. The gdb was locking itself when I was in an mxd that was accessing the feature classes inside it.

No other application was accessing it at the time i.e. ArcCatalog was closed. I would reboot the VM, open a blank map document, drag in the feature class I wanted, then try to edit that feature class, or create a new feature class in that gdb. No such luck.

I mentioned that this happened for about a week. After that, everything was fine again. I’m not sure why. Windows Updates have broken file gdbs before, but I don’t think it was that.

I think that IT changed some group policies and that affected my host where the data is stored, which in turn affected the VM. I know that sounds crazy, but I’m used to accepting crazy by now when it comes to Arc. It’s still not enough to make me go to QGIS though.

View list of ArcGIS Server Data Store Items

Recently I had a task where, for >10 ‘clients’, I had to:

1. Create a unique database user
2. Grant privileges to certain datasets in a SDE database
3. Save an mxd showing only what was relevant to that client
4. Publish a feature service for consumption in ArcGIS Online

While I will leave the specifics of that delightful script for another post, during debugging I came up against an error where the script would fail if the database connection had already been registered with ArcGIS Server.

	#
	# @date 23/03/2015
	# @author Cindy Williams
	#
	# Prints the list of databases registered with the
	# given ArcGIS Server, along with the connection
	# properties (excluding ENCRYPTED_PASSWORD).
	#
	# For use in the Python window in ArcCatalog.
	#
	import arcpy
	import os
	folder_user = os.environ['USERPROFILE']
	folder_arccatalog = "AppData\Roaming\ESRI\Desktop10.3\ArcCatalog"
	ags_name = "arcgis on dev01 (admin).ags"
	ags = os.path.join(folder_user, folder_arccatalog, ags_name)
	for dsi in arcpy.ListDataStoreItems(ags, "DATABASE"):
	print dsi[0] + "\n\t" + "\n\t".join(dsi[1].split(";")[1:])

view raw

viewAGSDataStoreDatabases.py

hosted with ❤ by GitHub

In lines 15 and 16, I get the location of the ArcCatalog connections folder for the current user. This is dependent on knowing which version of Arc is installed. I have found a slightly better way (written quite verbosely for clarity):

	import arcpy
	import os
	# Detailed description
	os_appdata = os.environ['APPDATA'] # Current user's APPDATA folder
	folder_esri = "ESRI" # ESRI folder name
	arc_prod = arcpy.GetInstallInfo()['ProductName'] # Get the installed product's name e.g. Desktop
	arc_ver = arcpy.GetInstallInfo()['Version'] # Get the installed product's version number
	arc_cat = "ArcCatalog" # ArcCatalog folder name
	print(os.path.join(os_appdata,
	folder_esri,
	arc_prod + arc_ver,
	arc_cat)

view raw

getArcCatalogConnectionsFolder.py

hosted with ❤ by GitHub

My dev environment

A while back, I made a big commotion about migrating my dev environment from Eclipse with PyDev to Visual Studio 2012 with PVTS. That did not last as long as I intended.

While there is nothing wrong with the setup, my scripts simply aren’t big enough or interwoven enough to warrant a gigantic IDE. Also, my complete failure in getting Git to play nice with it has led me to fully switch over to Notepad++.

I initially installed it to be able to quickly edit ArcGIS for Windows Mobile .amp files. I also used it when I was dabbling with HTML5. I then started noticing how convenient it was to open up my Python scripts in there, and how I could open a new file, scribble some pseudocode before leaving work, and having it persist between sessions.

I was hooked. I set it up with two views so I could view files side by side, played with a few themes before settling on Navajo, and  I was off.

Three problems though: I could not, for the life of me, run Python scripts from the editor. It was incredibly frustrating to have to do my edits, then open the file in Idle to run. Secondly, no arcpy intellisense (sad). Lastly, it does not recognise the .pyt extension of ArcGIS Python toolboxes.

Which is why I am now creating all my Python toolboxes and other scripts in PyScripter. I tweaked it to recognise .pyt files using Joel McCune’s technique, so now I can get some arcpy intellisense in my .pyt files now. I still miss my Notepad++ though, so I still find myself writing one-off scripts there, while leaving my bigger projects in PyScripter.

ArcGIS Diagrammer: Sad to see you go

I would often use ArcGIS Diagrammer for bigger projects, ones which would require documentation. I’ve never actually included any of my diagrams along with the project documentation, mostly because I was never asked for it.

Once I had worked with it long enough, I started to notice that it could be quite cumbersome and unintuitive at times. If something did not work as expected, I would first carry out the action manually in a file gdb, export the xml and look at it in Diagrammer to see the workflow.

Eventually, I just ended up going back to creating my templates in a template gdb. If I ever need a snapshot of my gdb in a Visio style, I will definitely open it up again, but it just doesn’t fit into my daily workflow.

What I am a bit irked about is the fact that ESRI never bothered to come up with a more official version of this product.

Coming up with a GIS Dataset naming convention

I have tried for several years to come up with/stick to a decent naming convention for GIS data. A colleague of mine came up with a 3 letter prefix scheme followed by an underscore which was pretty nifty:

• ftr – Feature class
• tbl – Table
• svw – Spatial view
• cvd – Coded value domain
• vew – Table view
• loc – Address locator
• gcs – Geocoding service
• svc – Service
• ras – Raster
• rvd – Range domain

This was done because when looking at the tables via SQL Server Management Studio, one cannot determine which tables are actually GIS tables/feature class attribute tables, and which are native SQL tables at a glance. This solved that problem, and conveniently grouped the different types of data at the same time.

However, a convention only works when everyone sticks to it. The use of plural vs singular is another bone of contention. I am of of the opinion that all feature classes should be named in the singular. This approach has stuck with me since mapwork classes in primary school (ah the 90s). Legends on maps always indicated “River”, “Road”, “Farm”, never “Rivers”, “Roads”, or “Farms”.

This is why when someone gives me shapefiles or whatever that have plural names, I immediately change it to singular. I WILL NOT TOLERATE THIS LUNACY. These naming issues fall under the greater problem of ever expanding formats for GIS data, along with the issue of “How on earth are we going to store all this stuff?!”

And no, saying “in the cloud” is not a solution. Someone still has to manage that data, and structure it in a way that makes sense and somehow anticipate further expansion of the system. That, however, is a story for another day.

ArcMap Woes Part 2: Assigning keyboard shortcuts

I came across this gem the other day. I was in the midst of a massive digitising session, ArcMap was being “slow” and my frustration levels were at a monthly high. I decided to spend a few minutes assigning keyboard shortcuts for the editing mode tools Edit Vertices and Continue Editing. What happened next just pushed my frustration levels higher – the shortcuts would not “assign”, no matter how many times I tried it, or what manner or key combination I used.

That GeoNet thread confirmed that there is a bug of some sort which prevents it from working properly. So I gave up, opened a blank mxd and carried on in there.