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.

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!