Dividing a polygon into a given number of equal areas with arcpy

I was recently searching the Internet trying to find any tool that would let me split an arbitrary polygon inside a geodatabase feature class into multiple polygons of equal area. ArcGIS does provide this functionality as a part of the parcel fabric functionality. Unfortunately, there is a lot of work involved in setting up the parcel fabric and there is a lot to learn before you will be able to divide your parcels. So I was looking for a simpler solution that would work directly with the geometry of a polygon.

However, I was not able to find any solution that would work and most helpful posts on the forums were pointing either at parcel fabrics or providing some ideas on the implementation of the workflow using multiple geoprocessing tools and some custom code. I found a nice ArcGIS custom script tool called Polygon Bisector which

computes a line that bisects, or divides in half, a polygon area along a line of constant latitude or longitude.

So, this would work great if you need to split a polygon creating a number that follows the exponent of two (2, 4, 8, 16, 32 and so forth). This is because after dividing a polygon into two polygons of equal area (now you have 2 polygons), you could divide each of them into two parts again (now you have 4 polygons), and so on. Since I want to be able to divide a polygon into an arbitrary number of areas, I had to write my own tool.

I have solved this problem this way. Say I want to have the polygon of area 1000 sq. m. divided into 5 equal areas:

  1. Get an extent of a polygon.
  2. Construct a polyline using the vertices of the polygon’s extent with a tiny shift of coordinates.
  3. Cut the polygon into two halves using this line.
  4. Find what is the area of the smallest polygon.
  5. If the area is smaller than the 200 sq. m. (that is, fifth part of the polygon), the shift the line again and re-run steps 2-4.
  6. If the area is 200 sq. m. or larger, than leave this part and keep working with the polygon that is left essentially running through the steps 2-5.
  7. When the original polygon has been successfully divided into equal areas, they are inserted into a new feature class along with the source polygon attributes.

The illustration of the cutting lines with the extent polygon is below.

extent_figure

This approach has several disadvantages, though.
First, if your polygon is very large and you want the parts of the polygon to have the same area with the minimal difference, the tool execution will take a lot of time as you will need to shift the cutting lines a few centimeters, cut the polygon into halves, and evaluate the result.
Second, you can only choose between the North-South or West-East direction of the cutting lines. You won’t be able to specify the angle yourself.

However, this tool works great for the target use I kept in mind when writing it. Using 0.5 meters as the step for moving the cutting line, the difference between the largest and the smallest sub-polygons was just around 1%. Running the same code for the same polygon with the shift value set to 0.05 meters (5 cm), I have observed the difference value to be around 0.1%.

The illustration of the polygon subdivision is below (West-East to the left, North-South to the right).

divide_polys_figure

The code is available at the GitHub Gist as usual. Run this code inside the Python window in ArcMap while having a single polygon selected.

 

 

Advertisements

Implementing Copy Parallel with ArcObjects

In ArcMap, in the Editor toolbar, you can find a command Copy Parallel which  will create a copy of a selected line feature at an offset distance you specify. This command is unique to the Editor toolbar and is not available as a geoprocessing tool. This makes it impossible to create parallel lines in batch (say, create lines to the both sides of the selected line feature starting with 10 meters up to 100 meters with the step of 20 meters).

A neat thing that exists in ArcMap is that you can assign a shortcut key to any command you can find in the Customize menu > Commands tab window. So, your workflow could look like this: you start the editing session, select a feature, use your shortcut key (e.g., Ctrl-Alt-P) and enter the offset value in the Copy Parallel window that appears. However, this implies that you would need to type an offset value for every copy to be created and you might also need to change other options in the window.

It’s important to understand that creating a parallel copy of a feature in the context of ArcGIS is not the same as creating a feature with a specified offset (because I thought it was). Compare the features produced using the Copy Parallel command (to the left) and using Python code (that just moves each vertex of a feature using an offset value) (to the right):

parallel

As you can see, the length of line segments can change when using the Copy Parallel command where as just moving vertices with offset preserves length of the segments.

You can call the Copy Parallel command using ArcObjects (just like any other command) but you won’t be able to supply any input parameters such as the offset value. So, it’s no different from assigning a shortcut key. The only solution left is to implement an own version of Copy Parallel that will mimic the built-in command. I have written an ArcGIS .NET add-in that can be called from a custom toolbar.

One can specify the intervals in a text box and then click the button.

parallelpara1

This is the code for the combo box item and for the button item:

 

 

Generate HTML report about Esri geodatabase using Python

Many GIS administrators and analysts maintain an inventory of all the data and metadata associated with a corporate Esri geodatabase. Most of the users have written short Python scripts that reads properties of geodatabase items and prints them to a text file or an Excel sheet.

A non-supported product called ArcGIS Diagrammer was able to generate HTML reports about your geodatabase reporting information about its domains, feature classes, tables and so much more. However, since version ArcGIS 10.3, the support for the tools was dropped, so if you want to use it now you need either to have ArcGIS Desktop 10.2 installed or have a bunch of .dll files from the ArcObjects .NET SDK 10.2 stored on your machine. This makes using ArcGIS Diagrammer cumbersome.

Anyways, ArcGIS Diagrammer wasn’t good enough for my needs because I have always wanted to work with an interactive representation (where I could sort columns, filter irrelevant information, and create printable views). Being frustrated with the need to hack own scripts for this, I’ve decided to write an own Python package that would report information about Esri geodatabase in an interactive HTML page.

Spending some evenings now and then, and that’s it – the package is in beta now, ready to install and use! Make sure to get it from its GitHub repository: registrant. This is how a sample report could look like; below is an animated example just to give you a feeling of what it would be like to work with the report.

report_sample

I would love to hear from users whether it’s useful and what other information and features such a report should have. You can submit an issue on the registrant‘s GitHub page or just leave a comment to this post.

 

Rounding double fields error in ArcMap

I have received an odd error message recently when trying to set up a label expression in ArcMap using Python parser.

I have a field of Double type and I want to show its value rounded in a label. Shouldn’t be more difficult than

int([AREA])

But I am getting an error message:

The expression contains an error.
Modify the expression and try again.

Error 0 on line 0.
Error running expression: esri__FindLabel(ESRIExpressionArg0)
Traceback (most recent call last):
File “<expression>”, line 1, in <module>
File “<string>”, line 2, in esri__FindLabel
ValueError: invalid literal for int() with base 10: ‘6380,614’

This is interesting. As it turns out, the double field’s string representation is used during the evaluation. On a virtual machine I’ve used, the locale that I’ve used had a comma , as the decimal symbol. So, in order to make it work I have to use the str.replace method:

int(float([AREA].replace(',', '.')))

Alternatively, if you want/can change the locale settings in Windows to use a dot . instead of a comma ,, then you will not need to run the replace method. It would be sufficient to run int(float([AREA]))

Small thing but hope this will save some time for other peer users!

Using SQL Server constraints and geodatabase domains

Attribute constraints

Many ArcGIS users use geodatabase domains to allow data editors to enter for certain attributes only certain values, either within a range or from a set of coded values. This functionality streamlines data editing and is very helpful indeed. However, having a geodatabase domain set for a field doesn’t actually prevent users from typing in value that is either outside of range or not in the list of coded values.

Entrance of illegal values can be done either programmatically using arcpy or SQL or by editing the attribute table of a feature class using Field Calculator. To be able to find out which features have illegal attribute values, you would need to select all of your features in the editing session and click Editor > Validate Features. This will select features with illegal values.

But what if you would like to let your users pick only certain values when editing the features and prohibit entering any illegal values? To do this, you could use database constraints such as foreign key constraint. In fact, I have already answered this exact question on GIS.SE: Restrict values to domain codes (beyond the attribute table).

In the end of the post, please look at the code snippet of what should be done in SQL. 

Now you can just use GP tool Table To Domain which will let you create a geodatabase domain from the dbo.CityType table as well as add the coded values into it. Then you can assign this domain to a field Type in the Cities feature class using the GP tool Assign Domain To Field.

Now user will get an error window in ArcMap (returned from SQL Server) every time they will try to enter illegal values into the field and save the feature. One thing to keep in mind when embracing this workflow is that you’d need to go to Editor toolbar > Options > Attributes tab and enable the option Display the attributes dialog before storing new features. This is necessary to do if you don’t specify any default value for this field that is stored within the dbo.CityType table. In this case, newly created features will have no value associated with the Type attribute and you won’t be able to digitize a feature on the map without getting the error message.

Spatial constraints

Another thing that may bug you is the geodatabase topology. It’s very handy when you have inherited a large suite of feature classes and you would like to enforce some integrity rules concerning the spatial relationships between features in those feature classes. However, if your data is stored in a multi-user geodatabase, then you could create own rules that would prohibit users from creating features that break those rules. Using ArcGIS geodatabase topology it is still possible to create a feature that would be considered invalid in terms of its relationship with another feature (say school point inside a lake polygon), however the only way to find this out is to validate topology on existing features.

Using SQL Server triggers, it is possible to specify the spatial rules and prevent creation of features that don’t follow these rules. Below is a simple example of a trigger that won’t let ArcMap users to digitize a point on the map to create a point feature that is located outside of the boundaries of the California state.

Network Analyst and ArcPy: finding k-alternate path

Using ArcGIS Network Analyst extension, it is possible to find out the best route between a pair of points (best is defined in terms of the impedance — it can the shortest or the fastest route, for instance). However, you might like having multiple alternatives to the best route which one usually refers to as the K shortest paths.

The K shortest path routing algorithm is

an extension algorithm of the shortest path routing algorithm in a given network.

When solving a route between a pair of points, Network Analyst finds only a single route that is considered to be the best. That is, one cannot get a number of routes while solving a Route layer. However, this is something many of you have seen when routing using Google Maps or HERE, for instance.

I got curious and have started searching for this algorithm implementation using ArcGIS Network Analyst. I’ve found a very useful Esri forums post by a former Esri Network Analyst product engineer Michael Rice who wrote that

At the ArcGIS 10 release, you can approximate a k-shortest paths solver using the enhanced barrier features. Additionally, for line and polygon barriers, instead of just restricting the parts of the network that are covered by their geometry, you can opt to simply scale the costs of the parts of the network covered by the geometry.

Therefore, to approximate k-shortest paths, you can do the following:

1. Solve the route to get the best path
2. Take the current best path and load it as a polyline barrier with a scale factor of x
3. Repeat

Note that by scaling the polyline barrier instead of restricting it, this will allow
the next path to reuse some of the same edges along the previous path if necessary, but at a higher cost than in previous paths. Also note that when multiple line barriers apply to the same edge, the final scale factor for the cost of that edge becomes the product of their individual scale factors. For example, if there are two line barriers applied to some edge with scale factors x and y, respectively, the scaled cost of traversing that edge will be x*y*originalEdgeCost.

This means that the more an edge is reused across different paths, the less likely it will be to be used in any subsequent paths as the process discussed above continues to iterate.

This approach will technically not guarantee a true k-shortest path solve, so we must
be careful not to incorrectly label it as such. As I said previously, it is merely an
approximation for this concept. A more appropriate name might be to call this a
k-alternate path approach (as opposed to shortest). The alternate paths you get will
ultimately depend on the scale factors you use.

I’ve written a Python add-in that given a number of alternative routes to find, will generate a new feature layer with a number of routes solved. I have used the workflow Michael described loading a solved route as a polyline barrier for all subsequent solve operations. You can use this code when you want to generate multiple alternate routes. This is how it looks in ArcMap.

2017-04-27 15_12_00-Routing.mxd - ArcMap

The Python code to be used to make the add-in:

ArcGIS Pro vs ArcGIS Desktop: geoprocessing tools

Both ArcGIS Desktop and ArcGIS Pro have a suite of toolboxes with hundreds of geoprocessing tools available. However, there are still quite a few tools that are available only in ArcMap and have not been ported to Pro. You can find the list of those tools at the Help page Tools that are not available in ArcGIS Pro.

However, what not so many people may know is that there are actually quite a few tools that are available exclusively in ArcGIS Pro. They have never existed in ArcMap in other words. They may or may not become a core part of the ArcGIS Desktop ArcToolbox. I have seen that many of the tools added to ArcGIS Pro 1.3 have been added to ArcGIS Desktop 10.5. A good part of them is specific to Pro, however, and would not make sense in ArcMap, such as Consolidate Project or Package Project tools. However, there are some tools that are applicable to both environments as they operate on geodatabase, for instance, Alter Domain, or with a network dataset, for instance, Create Network Dataset From Template.

To be able to find out what tools are available in ArcGIS Pro and not available in ArcGIS Desktop, I’ve written a tiny Python script that reads all the toolboxes available in the both installations and then compares them. You would need to run the code for Pro and Desktop Python installations, of course.

As I need to call Python installed in Pro and Desktop, I prefer to use a .bat file for this:

"C:\Python27\ArcGIS10.4\python.exe" "C:\GIS\Temp\ToolLister.py" "Desktop"
"C:\Program Files\ArcGIS\Pro\bin\Python\Scripts\propy.bat" "C:\GIS\Temp\ToolLister.py" "Pro"

After calling this .bat file, you will get two .json files each containing the tools available in these two environments.

Please note that I have run this code for ArcGIS Desktop 10.4 and ArcGIS Pro 1.4.

The code for your reference (ToolLister.py):

The list of tools present in Pro that are not in Desktop

ArcGIS Pro does have some new toolsets we never had in Desktop.

Standard Feature Analysis

The Standard Feature Analysis toolbox contains tools for performing spatial analysis on feature data in your portal.

Raster Analysis

The Raster Analysis toolbox contains a set of powerful tools for performing raster analysis on data in your portal. By distributing the processing between multiple server nodes, you can efficiently process large datasets in less time than ever before.

GeoAnalytics Tools

The GeoAnalytics Tools toolbox contains a set of powerful tools for performing spatial analysis on big data. GeoAnalytics Tools are powered by your ArcGIS GeoAnalytics Server. ArcGIS GeoAnalytics Server distributes the analysis between multiple server nodes. By using distributed processing, you can process large datasets in less time.

And here are the new tools in existing toolboxes:

Geostatistical Analyst Tools
  • EBKRegressionPrediction_ga
  • GALayerToRasters_ga
Editing Tools
  • AlignFeatures_edit (since Desktop 10.5 only)
  • CalculateTransformationErrors_edit (since Desktop 10.5 only)
Analysis Tools
  • EnrichLayer_analysis
  • PairwiseIntersect_analysis
  • SummarizeWithin_analysis
  • PairwiseDissolve_analysis
  • SummarizeNearby_analysis
  • PairwiseBuffer_analysis
  • SplitByAttributes_analysis
  • GraphicBuffer_analysis
Spatial Analyst Tools
  • LocateRegions_sa (since Desktop 10.5 only)
  • RemoveRasterSegmentTilingArtifacts_sa (since Desktop 10.5 only)
3D Analyst Tools
  • ClassifyLasBuilding_3d (since Desktop 10.5 only)
  • TileLas_3d (since Desktop 10.5 only)
Conversion Tools
  • AddRasterToGeoPackage_conversion (since Desktop 10.5 only)
Data Management Tools (some of the tools may be available since Desktop 10.5)
  • AlterDomain_management
  • AddRelate_management
  • RemoveRelate_management
  • CreateMobileMapPackage_management
  • CreateVectorTileIndex_management
  • CreateVectorTilePackage_management
  • MosaicDatasetToMobileMosaicDataset_management
  • CreateSceneLayerPackage_management
  • EnableCOGO_management
  • DisableCOGO_management
  • BuildStereoModel_management
  • InterpolateFromPointCloud_management
  • ComputeCameraModel_management
  • GeneratePointCloud_management
  • ComputeMosaicCandidates_management
Space Time Pattern Mining Tools
  • LocalOutlierAnalysis_stpm
Network Analyst Tools
  • CreateNetworkDatasetFromTemplate_na
  • CreateTemplateFromNetworkDataset_na

Hope you will find this list useful. If your main workflows are still ArcGIS Desktop based (which I believe they are for most of us), keep in mind that you can call ArcGIS Pro Python 3.5 installation to take advantage of those new geoprocessing tools if you have both Pro and Desktop installed. Take a look at this Help page Python in ArcGIS Pro. To learn more about calling Windows applications from Python, take a look at the built-in subprocess module.