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:

 

 

Installing GDAL with Anaconda on Windows

It has been historically fairly hard to install GDAL along with all its C dependencies and make it play nicely with existing Python libraries leaving alone having multiple environments with multiple versions.

I have blogged a year ago on how to set up a nice open-source GIS sandbox using Anaconda. Nowadays, it got even easier. When installing Anaconda 4, an application called Anaconda Navigator will be installed. From the web site linked:

Anaconda Navigator is a desktop graphical user interface (GUI) that is included in Anaconda® and allows you to launch applications and easily manage conda packages, environments and channels without using command-line commands. You can configure Navigator to search for packages on Anaconda Cloud or in a local Anaconda Repository. It is available for Windows, macOS and Linux.

Using the Anaconda Navigator, you can install GDAL without even starting a terminal!

Here are the steps you need to follow:

  1. Download and install Anaconda 4.x, 64-bit, Python 2.7, Windows version.
  2. Make sure to install Microsoft Visual C++ 2008 Service Pack 1 Redistributable Package. If you don’t do that, you may get weird errors when trying to import gdal in your Python code.
  3. Install gdal 2.1.0 using the Navigator application.
  4. Optionally, add a system environment for specifying the gdal folder in your Anaconda installation folder –GDAL_DATA =  C:\Users\dev\Anaconda2\Library\share\gdal – if you would need to work with spatial references, for instance, use pyproj to re-project coordinates.
  5. Start a Python command prompt and run import gdal and import ogr. If both have been successfully imported, you are good to go!

It is very easy to manage your packages using Navigator. Make sure to add some other Conda channels such as conda-forge because you won’t be able to find all geospatial packages you might need in the default channel. Thereafter, you can easily add pyprojand other geospatial Python packages such as pysal, descartes, basemap and cartopy.

Read file geodatabase domains with GDAL and use XML schema

As I keep working on the registrant Python package, I have decided to add support for generating report about file geodatabase for a case when you don’t have any ArcGIS software installed on a machine yet would like to use the package.

To be able to read an Esri file geodatabase using GDAL, you can either use OpenFileGDB driver or FileGDB driver which relies on FileGDB API SDK. I thought that it would be nice to avoid being dependent on a third-party library, so I am using plain GDAL which makes it possible to interrogate a file geodatabase without having any ArcGIS software installed.

You can look into the source code to see how information about tables and feature classes as well as their fields can be pulled. In this post, I will just show you how you can read file geodatabase domains. This is a bit special because there are no built-in methods for reading domains. There is a GDAL enhancement ticket that targets this, but for now the only workaround I found is to run a SQL query against a metadata table, get an XML string back and then parse it. Esri has published a technical paper XML Schema of the Geodatabase which definitely helps to navigate the XML representation.

Using domains can be handy in a situation when you would like to report not the actual data stored in a table and accessible through OGR (that is, codes), but rather the human-readable representation (that is, values). For doing this, you would need to grab an XML representation of a particular feature class and see what domains are used for a specific field.

The Python code for reading domains and then finding out what fields have domains assigned is below:

 

 

Geolocation and mapping with Elasticsearch

I have played around with Elasticsearch for a while and it has been my first time I was working with a NoSQL database. It’s been a lot of fun to load some data and see what querying techniques are available. I was particularly interested in seeing what kind of support for geospatial operations does the elasticsearch provides. As it turned out, there is very good support for GeoJSON data structures for storage and visualization (both as point data and as areas). You can even run the spatial queries (e.g. find all points within a specified polygon) and draw the results on a tiled basemap as an operational layer (using Kibana) along with using regular SQL-like attribute data filtering.

You can load your geospatial datasets into elasticsearch NoSQL database and then visualize it on a basemap as well as run some spatial querying, draw heat maps and choropleth maps using own regions (polygons that can be enriched with the point data). This is just a few maps I’ve created really quickly loading the cities and states geodatabase feature classes into the NoSQL database and then using Kibana web interface to author visualizations:

It all seems to be very powerful yet fairly easy to configure. Loading the data into the database was not very straightforward, though, because as I understood, you cannot feed the source GeoJSON files into the elasticsearch when using the _bulk api to load documents into an index. You can read more about GeoJSON support in elastic in the docs page Geo-Shape datatype.

I’ve written a tiny Python script that will convert a geodatabase feature class into a GeoJSON and then construct a new .json file with the proper format that can be loaded into the elasticsearch database.

To play with the elasticsearch API interface, I’ve used Postman. It’s a very handy application that will let you construct all kinds of HTTP requests and save them to reuse later on. Using Postman, I was able to submit the PUT request to load documents from the .json file I’ve created by running the Python script using the binary body and browsing to the source data file.

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.

 

Convert an Esri geodatabase feature class to GeoJSON

I bet many of you using Esri geodatabase datasets needed at some point to convert your data to some other formats for interoperability. GeoJSON is one of those formats that can be understood by many other systems and APIs.

If you have ArcGIS Desktop 10.5+, you can use Features To JSON geoprocessing tool as in version 10.5, Esri has added optional geoJSON parameter. Please mind that this tool didn’t have this parameter in previous versions.

If you are still on ArcGIS Desktop 10.4 or an earlier version, you could use the private interface of arcpy.Geometry objects that can expose some of the information about the geometries in the form of geoJSON and then do the rest on your own using json module:

g = arcpy.PointGeometry(arcpy.Point(45, 45), arcpy.SpatialReference(4326))
g.__geo_interface__
{'coordinates': (45.0, 45.0), 'type': 'Point'}

If you have ArcGIS Pro installed, then you can use the Features To JSON geoprocessing tool which has always had this optional geoJSON parameter which will do the trick.

If you need to do a conversion on a machine without ArcGIS Desktop, you could use ogr2ogr command (which will convert the Counties shapefile into a .json file in geoJSON format):

.\ogr2ogr.exe -f “GeoJSON” “C:\GIS\output.json” “C:\GIS\Temp\Counties.shp” -select NAME,POP2000

You could also use some other conversion tools such as ogr Python module:

What I really like about geoJSON is that it’s possible to draw its features online very easily such as using http://geojson.io (you can even edit, style and share your geoJSON using this service!) or on GitHub (they have pretty good support for geoJSON). On GutHub, you can even show the geoJSON features on a basemap directly in the source code or in a gist. Look at this example!