Module 2D, GeoDMS and Python

learning objective: exchanging tabular data between GeoDMS and Python (pandas), and knowing which tool to use for which part of your model

introduction

Many GeoDMS modellers also work in Python. The two are not competitors: they are good at different things, and this module shows how to make them work together.

The core of the model belongs in GeoDMS: every data item has explicit units that the engine checks, calculation is lazy (only what you request is calculated), and every result is traceable back through its expression and suppliers. Python belongs at the edges: estimating parameters before the model runs, machine learning, bespoke statistics, and making figures and reports from the results afterwards. The main wiki page Why GeoDMS explains this division of labour in detail and contains a translation table from pandas/numpy concepts to GeoDMS concepts (DataFrame index = domain unit, column = attribute, groupby-sum = aggregation with a relation, and so on). If you know pandas, that table is the fastest way to map your existing knowledge onto what you learned in Module 1.

There are three proven architectures for combining the two, described on GeoDMS through Python:

  1. Python drives: run the GeoDMS as a separate process from Python.
  2. GeoDMS drives: run a Python script as a separate process from within the GeoDMS.
  3. In-process: embed the GeoDMS engine in Python through the geodms module.

The hands-on part of this module covers the data exchange that architectures 1 and 2 are built on: writing and reading tabular data as parquet files. At the end we briefly look at the other routes.

why parquet?

Apache Parquet is a columnar, binary file format for tabular data. Unlike CSV, it stores the value type of each column in the file itself and compresses the data. That makes it the recommended format for exchanging tables with Python: an int32 column written by the GeoDMS arrives in pandas as an integer column, and vice versa, with no conversion code on either side.

The GeoDMS reads and writes parquet through the GDAL library, with the same gdal.vect and gdalwrite.vect StorageManagers you used for CSV and shapefiles in Modules 2A and 2B. The .parquet file extension is enough for the GeoDMS to select the GDAL Parquet driver automatically. See the main wiki page Parquet for the full reference.

One limitation to keep in mind: attributes of all value types are written, except the point group. So parquet exchange is for tables; keep the geometry on the GeoDMS side and relate the Python results back to the spatial domain, as we do below.

On the Python side you need pandas with a parquet engine, usually pyarrow:

pip install pandas pyarrow

from GeoDMS to pandas

We continue with the data_sources project from Modules 2A and 2B, where the SourceData container reads gemeente.csv (Dutch municipalities) into the domain unit gemeente. Recall from Module 2A that all CSV columns are read as string by default, so we convert the number of inhabitants (AANT_INW) to a numeric type in the export.

Writing a parquet file follows exactly the export pattern from Module 2B, only with a .parquet extension:

container Export
{
   unit<uint32> gemeente_export := SourceData/gemeente
   ,  StorageName     = "%LocalDataProjDir%/python/gemeente.parquet"
   ,  StorageType     = "gdalwrite.vect"
   ,  StorageReadOnly = "False"
   {
      attribute<string> GM_CODE  := SourceData/gemeente/GM_CODE;
      attribute<string> GM_NAAM  := SourceData/gemeente/GM_NAAM;
      attribute<int32>  AANT_INW := SourceData/gemeente/AANT_INW[int32];
   }
}

Trigger the write by double-clicking gemeente_export in the GeoDMS GUI, or by right-clicking it in the TreeView and choosing Update Subtree (Ctrl+T).

Now read the file in Python. Three lines of pandas are enough to load the table and list the five largest municipalities:

import pandas as pd

df = pd.read_parquet("C:/LocalData/data_sources/python/gemeente.parquet")
print(df.nlargest(5, "AANT_INW")[["GM_NAAM", "AANT_INW"]])

%LocalDataProjDir% expands to a project-specific folder under your LocalData folder (by default C:/LocalData/<project folder name>). Check where the file appeared on your machine and adjust the path in the Python script accordingly.

Note what did not happen: no type guessing, no separate metadata file, no conversion code. AANT_INW arrives in pandas as an integer column because the value type travelled inside the parquet file.

from pandas back to GeoDMS

Now the other direction. Suppose Python computes something new, here simply a rank by number of inhabitants, and writes it back:

df["INW_RANG"] = df["AANT_INW"].rank(ascending=False).astype("int32")
df[["GM_CODE", "INW_RANG"]].to_parquet(
    "C:/LocalData/data_sources/python/gemeente_python.parquet", index=False)

Use index=False, otherwise pandas writes its row index as an extra column.

Reading the file into the GeoDMS creates a new domain unit. If your script did not sort, filter or add rows, the rows still correspond one-to-one and in the same order to the existing gemeente domain, and you can tie the read domain directly to it with a calculation rule:

unit<uint32> gemeente_python := SourceData/gemeente
,  StorageName     = "%LocalDataProjDir%/python/gemeente_python.parquet"
,  StorageType     = "gdal.vect"
,  StorageReadOnly = "True";

Note that the unit body can stay empty: all columns in the parquet file become available as attributes automatically, with their value types taken from the file.

Relying on row order is fragile, though. The moment your Python script filters rows or sorts the DataFrame, the one-to-one correspondence silently breaks. The safe route is to relate both domains on an external key, here the municipality code, with rlookup and the -> shorthand of lookup:

unit<uint32> gemeente_python
:  StorageName     = "%LocalDataProjDir%/python/gemeente_python.parquet"
,  StorageType     = "gdal.vect"
,  StorageReadOnly = "True";

attribute<gemeente_python> py_rel   (SourceData/gemeente) :=
   rlookup(SourceData/gemeente/GM_CODE, gemeente_python/GM_CODE);
attribute<int32>           INW_RANG (SourceData/gemeente) := py_rel -> INW_RANG;

Municipalities whose code is not found in the parquet file get a null relation, and therefore a null INW_RANG. That is exactly what you want: missing data stays visibly missing instead of being silently misaligned. The Parquet page describes a third option, copying values one-to-one with union_data, which raises an error when the row counts differ.

In this module you ran the Python script by hand between the write and the read. The GeoDMS can also start the script itself, automatically, at the moment its output is requested, by weaving the ExitCode of the exec_ec function into the StorageName of the output. That page contains a worked end-to-end example of the write, run, read chain.

the other two integration routes

Python drives: GeoDmsRun. When Python orchestrates the whole workflow (for example a calibration loop, or a batch of scenario runs), it can start the command line tool GeoDMSRun (GeoDmsRun.exe) with subprocess to update one or more items of a configuration. The GeoDMS writes its results to the storages configured in the model, for example the parquet exports you built above, and Python reads them back. This requires no GeoDMS-specific Python code and keeps both processes fully isolated. See architecture 1 on GeoDMS through Python for a code example.

In-process: the geodms module. The GeoDMS also ships a geodms Python module (PyDms.pyd) that embeds the engine inside the Python process. Python can then load a configuration, navigate the item tree, change expressions, update items and read values back, without starting a separate process and without exchanging files. The bindings are a work in progress and tied to specific Python versions; see Python bindings for the version notes and the full API reference.

try it yourself!

  • Download the project here if you have not done so already, and unzip it to C:/prj/GeoDMSAcademy.
  • Open exercise.dms (in the GeoDMS_Academy/data_sources/cfg subfolder) in your text editor, with your SourceData configuration from Module 2A in place (or start from result.dms in the same folder).

Now build a full roundtrip:

  1. Add an Export container that writes GM_CODE, GM_NAAM and AANT_INW (converted to int32) of the gemeente domain to a parquet file, as shown above. Trigger the export in the GeoDMS GUI and verify the file appears in your LocalData folder.
  2. Write a small Python script that reads the parquet file with pandas, removes the rows with missing data (in this dataset missing values are coded as -99999999), computes for each remaining municipality its share of the total number of inhabitants, and writes GM_CODE plus the new share column to a new parquet file.
  3. Read the result back into the GeoDMS. Because your script filtered rows, the one-to-one route is not safe anymore: use the rlookup route on GM_CODE to bring the share onto the gemeente domain. Open the attribute as a table and check that the filtered-out municipalities show null values.

As a bonus, think about which parts of this roundtrip you would automate with exec_ec if the Python step were a permanent part of your model, and read the worked example on the exec_ec page.

Go to previous module: Module 2C, GeoDMS own data formats

Go to next module: Module 3, Meta scripting ‐ templates and for_each