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:

 

 

Advertisements

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!

 

 

 

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 analysis with Esri routing services using Python

Good news for everyone who loves Python and routing! Now you can submit the routing requests using your input stops with the help of Esri Directions and Routing Services available to all organizations with the ArcGIS Online subscription and the ArcGIS Python API.

Maybe you use the ArcGIS API for JavaScript based web application or a Web App Builder application with the Directions widget or the ArcGIS Online Map Viewer. Some other users have built own Python (or other language based) scripts accessing the Esri routing services via the ArcGIS REST API in an own script or an application. I have been doing this using the ArcREST Python package.

Now you can use the ArcGIS API for Python to submit synchronous or asynchronous calls to the various routing services. Please refer to this ArcGIS REST API Help page to learn more what solvers are supported and how each of them works.

With the new Python API, working with the routing services became really easy. Many sample notebooks showcasing the routing services just got published to the ArcGIS API for Python documentation page.

Keep in mind that some of the routing related tools are also available as a part of Spatial Analysis Services. The notebook for these tools is available here.

Synchronous solvers:

Asynchronous solvers: