Printing pretty tables with Python in ArcGIS

This post would of interest to ArcGIS users authoring custom Python script tools who need to print out tables in the tool dialog box. You would also benefit from the following information if you need to print out some information in the Python window of ArcMap doing some ad hoc data exploration.

Fairly often your only way to communicate the results of the tool execution is to print out a table that the user could look at. It is possible to create an Excel file using a Python package such as xlsxwriter or by exporting an existing data structure such as a pandas data frame into an Excel or .csv file which user could open. Keep in mind that it is possible to start Excel with the file open using the os.system command:

os.system('start excel.exe {0}'.format(excel_path))

However, if you only need to print out some simple information into a table format within the dialog box of the running tool, you could construct such a table using built-in Python. This is particularly helpful in those cases where you cannot guarantee that the end user will have the 3rd party Python packages installed or where the output table is really small and it is not supposed to be analyzed or processed further.

However, as soon as you would try to build something flexible with the varying column width or when you don’t know beforehand what output columns and what data the table will be printed with, it gets very tedious. You need to manipulate multiple strings and tuples making sure everything draws properly.

In these cases, it is so much nicer to be able to take advantage of the external Python packages where all these concerns have been already taken care of. I have been using the tabulate, but there are a few others such as PrettyTable and texttable both of which will generate a formatted text table using ASCII characters.

To give you a sense of the tabulate package, look at the code necessary to produce a nice table using the ugly formatted strings (the first part) and using the tabulate package (the second part):

The output of the table produced using the built-in modules only:

builtin

The output of the table produced using the tabulate module:

tabulate

 

 

Advertisements

Multiple Ring Buffer with PostGIS and SQL Server

Recently I needed to generate multiple ring buffers around some point features. This can be done using a dozen of tools – Multiple Ring Buffer geoprocessing tool in ArcGIS, using arcpy to generate multiple buffer polygons and merging them into a single feature class using the buffer() method of arcpy.Geometry() object, or by using open source GIS tools such as QGIS. This is also possible to achieve using relational database that has support for the spatial functions. In this post, I would like to show you how this can be done using the ST_Buffer spatial function in PostGIS and SQL Server.

In order to generate multiple buffer distance values (for instance, from 100 to 500 with the step of 100) in SQL Server, I would probably need use CTE or just create a plain in-memory table using declare; in other words, this is what it takes to run range(100, 501, 100) in Python.

In the gist below, there are two ways to generate multiple buffers – using the plain table and the CTE.

Generating a sequence of distances in Postgres is a lot easier thanks to the presence of the generate_series function which provides the same syntax as range in Python.

Visualizing computational geometry concepts using JTS TestBuilder

In this post, I would like to let you know about an excellent piece of software, Java Topology Suite (JTS).

JTS is an open source library of spatial predicates and functions for processing geometries. It provides a complete, consistent, and robust implementation of fundamental algorithms for processing linear geometry on the 2-dimensional Cartesian plane.

A funny thing about it is that JTS

is used by most java based Open Source geospatial applications, and GEOS, which is a C++ port of JTS, is used by most C based applications.

So, all the downstream projects using GEOS such as various Python wrappers around GEOS such as shapely and even the PostgreSQL extension, PostGIS, all of them really work against the JTS using the GEOS as the interface for communication. So the JTS is a very, very powerful Java library.

If you are not a Java developer, though, this might be of little interest to you. However, there is another little application, called JTS TestBuilder, which provides a GUI for geometry exploration and is an interface into the JTS API. It is not so famous as other pieces of open source GIS stack, such as QGIS or GRASS, though. Also its documentation is outdated and scarce, so you would need to find out how to use the application on your own.

Nevertheless, it is an indispensable tool for anyone who spends a fair amount of time working with computational geometry or spatial data processing applications. It would also serve as a great visualization tool for GIS instructors who need to visually explain how GIS algorithms operate. I have used it to show how Convex Hull is created from a set of points, for instance. One obvious advantage of JTS TestBuilder is that you do not need to run any heavy GIS applications and the “geometry modification – running analysis – seeing the result” cycle is really short.

Here I’ve loaded cities of California along with the state boundary and created a convex hull for the boundary geometry.

2018-03-14 17_41_18-JTS TestBuilder

Having said that, you can work in the following manner:

  • Use your favorite GIS database management tool to get WKT of a geometry you would like to inspect or analyze.
  • Use the JTS TestBuilder to draw the features.
  • Run JTS Geometry Functions constructing new geometries or answering spatial questions.
  • Load the results of the analysis back into your GIS (either for ad hoc exploration or for storage).

The code to read the features into WKT and write back from WKT using arcpy:

 

JTS TestBuilder can also help you to learn something new with regard to GIS theory. If you think that you are a well seasoned GIS professional who can amaze others by mentioning a few cool names like Voronoi or Thiessen, I encourage you to explore the geometry functions JTS TestBuilder provides. I am pretty sure just a few of you have heard of:

  • Koch snowflake which are used a lot in space-filling as well as cartographic simplifaction algorithms.
  • Seirpinski carpet which is not used extensively in GIS yet, but there are some emerging applications regarding urban pattern analysis.

If you would like to take advantage of the computational geometry algorithms implemented in JTS, there are also ports to .NET and JavaScript.

Another very similar application that is particularly popular among math teachers is GeoGebra. I have been using it for a while, too, but it lacks export of result geometries into WKT which can be put into a geospatial database or drawn directly in a desktop GIS application such as ArcMap or QGIS. You can try GeoGebra online or by installing a desktop application. It is also available as an app for iOS, Android, and Windows Phone.

Desktop PyQt application for executing SQL queries against Esri file geodatabase

As I was always interested in building GUIs for some GIS operations, I thought that exploring PyQt deeper would be fun. A project has started as an experimental playground to see what functionality PyQt provides. As I spent more time working with PyQt, I have started wondering what it would take to build a useful desktop application.

Because I often find myself in need of querying a file geodatabase’s datasets, I have decided to build a GUI-based SQL editor that would let me execute SQL queries against a table or a feature class and draw the result set in a table form for further visual inspection. I have thought that other GIS users and developers may find this application useful and I therefore have decided to start a GDB: GitHub repository to let others take advantage of my work. Here it comes, check it out!

GDBee_sample

GDBee is a PyQt5 desktop application which you can use to write and execute SQL queries against tables and feature classes stored inside an Esri file geodatabase. It provides a tabbed interface which lets you connect to multiple geodatabases within a single session. It has a rich code editor featuring auto-completion (with suggestions), syntax highlight, and result set export.

If you are a QGIS Desktop user, you are already able to execute SQL against file geodatabases using QGIS DBManager plugin, but GDBee has some extra features that the DBManager is missing (for instance, you do not need to add your datasets as layers first and you can choose to copy individual cells instead of the whole row) from the result table.

Because Python is so widely used in the GIS community, I thought it would make sense to take advantage of Python bindings of GDAL (via GEOS) to be able to connect to a file geodatabase and execute SQL queries. Working with a file geodatabase via GEOS makes it possible to take advantage of SQL spatial functions that are otherwise inaccessible to an ArcGIS user!

The application provides multiple features:

  • Working with multiple geodatabases using multiple tabs (single geodatabase connection per tab)
  • Exporting result sets into various formats (WKT strings to paste into QGIS using QuickWKT plugin, arcpy code to paste into ArcMap Python window, pandas data frame via .csv file (which can be taken into geopandas), and Markdown table via .md file or plain text)
  • Executing SQL query with respect to the user selection (only selected text is executed)
  • Loading/saving SQL queries from and to text files on disk
  • Convenient keyboard shortcuts for query execution (F5 and Ctrl-Enter) and tab interaction (Ctrl-N and Ctrl-W for opening and closing tabs)
  • Copying data from the result set table (either individual cell values or row(s) with the headers preserved) – ready to paste properly into an Excel sheet
  • Choosing whether you want to have geometry column in the result set as WKT

You can look at its GitHub repository: GDBee: GitHub repo. You may find this PyQt desktop application useful if:

  • You would like to be able to interrogate your file geodatabase datasets using SQL (instead of Python-based interface such as Esri arcpy or open-source ogr)
  • You are an ArcGIS user that does not want to have QGIS Desktop installed just to be able to execute SQL against a file geodatabase
  • You use SQL on a daily basis working with spatial databases (such as PostgreSQL or Microsoft SQL Server) and want to be able to execute ad hoc SQL queries against file geodatabase datasets without loading them into a proper DBMS database
  • You already have a lot SQL code targeting tables stored in a DBMS spatial database and you would like to be able to reuse this code when targeting a file geodatabase

Do you think there is some other functionality that should be added? Please let me know by submitting an issue in the repository.

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.