This is a complete re-write using HighFive C++ bindings for HDF5. Issues on OSX are hopefully now resolved.

For more information see the list of changes.

For those building the in-tree GDAL driver, you will need the changes in GDAL PR 14596. Please open an issue for any problems that are found.

One last piece of TuiView functionality that has not yet been covered in this blog series is the ability to add “surrogate” colour tables to a file and use them to display images within TuiView.

Normally, a thematic file is displayed using the colour table within the Raster Attribute Table. However, in situations where RAT columns have different classifications which require the image to be shown in the colour table for that classification it is handy to be able to add surrogate colour tables and view the image with those.

The tuiviewwritetable utility

Installed with TuiView, the tuiviewwritetable command line utility allows surrogate colour tables to be added, removed and queried:

$ tuiviewwritetable -h
usage: tuiviewwritetable [-h] [-s SOURCE] [-n NAME] [-d DEST] [-p PRINTCT] [-r REMOVE]

options:
  -h, --help            show this help message and exit
  -s SOURCE, --source SOURCE
                        File to read color table from
  -n NAME, --name NAME  name to save the color table under
  -d DEST, --dest DEST  destination file to write color table into
  -p PRINTCT, --print PRINTCT
                        print out available color tables
  -r REMOVE, --remove REMOVE
                        remove table from specified file (must specify --name also)

Adding a surrogate colour table

Because surrogate colour tables are written to the image metadata, we recommend that only smaller (< 100,000 rows) tables are written using this method. Image metadata is much less efficient at storing the table information than the RAT.

When adding a surrogate colour table from another file, specify the --source, --dest and --name parameters. If more than one colour table is added using this method, TuiView will ask you which colour table you would like to use.

The following command will add a new surrogate colour table named “biomass” into myimage.kea from the image biomassimage.tif:

$ tuiviewwritetable --dest myimage.kea --source biomassimage.tif --name biomass

Querying surrogate colour tables

If you wish to see which colour tables are available use the --print option:

$ tuiviewwritetable --print myimage.kea
Name	    Size
------------------------
biomass	    65537
water       65537

Deleting surrogate colour tables

To remove a surrogate colour table, use the --remove flag in combination with the --name option:

$ tuiviewwritetable --remove myimage.kea --name biomass2

Accessing surrogate colour tables within TuiView

If you load a thematic image with TuiView, it will default to using the colour table within the Raster Attribute Table. However, if you open the Query Window and right click on a float or integer column you will see there is an option to “Set column as Colour Table Lookup”:

If you select this option, the behaviour will be either:

1. If you have only one surrogate colour table in the current image this column will be directly used as a lookup into the surrogate colour table.
2. If you have more than one surrogate colour table you will be presented with a dialog allowing you to choose which surrogate colour table to use. This column will then be used as a lookup into the chosen surrogate colour table.

The image will now be displayed using the column as lookup into the surrogate colour table:

Note that a small icon () is now drawn against next to the column header that is being used as a lookup.

To remove the colour table lookup, select the same menu option again and it will be removed and the original colour table shown.

Conclusion

TuiView surrogate colour tables allow you to view your imagery using multiple colour tables for different classifications. This can be a handy feature in situations where you have many different classifications stored in your Raster Attribute Table.

This is mainly a bugfix release. There have been some fixes to the way TuiView decides an image is thematic, vector labels are now always within their polygon, and new viewers always now show the current geolinked extent. There has been a fix for the tuiview process not exiting on Wayland, however other Qt/Wayland problems remain such as dockable windows being unable to be moved. Until this is fixed in Qt, using X with TuiView is recommended under Linux.

For more information see the list of changes.

Most modern computers come with a memory cache. This holds a copy of memory chunks most recently used and can be up to 100 times faster than normal memory. There are in fact often multiple levels of cache but in this discussion we just assume there is just one.

Each time you access memory that is not in cache (known as a cache miss), the CPU must go to main memory to retrieve the requested item. As it does this, it makes a prediction that the items immediately following will be the next ones requested, and loads a block of those into the on-chip cache at the same time (with very little extra overhead). If the prediction turns out to be correct, this means that those subsequent requests are met from the (much faster) on-chip cache, rather than requiring more requests to main memory.

In order to gain the most benefit from this, we should generally try to ensure that our processing matches the predictions the cache system is making, and so we should try to process in the order in which data is being held in main memory.

Directly monitoring this aspect of performance tends to be quite difficult. Confounding the analysis is the fact that operating systems tend to show that the CPU is ‘busy’ while waiting for data to come from main memory. The only way to determine whether a particular implementation is fast or not is to run it and time how long it takes. Using how ‘busy’ the CPU is might not give you an accurate picture.

A note about multi threading: Run multiple threads at once 
generally makes it harder for the cache to run efficiently 
since the memory requests become less predicable. Software also 
tends to have to add more locking and checking when multiple
threads are enabled and this overhead can affect the speed
of your code.

The below assumes there is no other competing processing running on your hardware, results will vary depending on what other things are going on in your system.

This simple 1D case

Assume you are looping over a 1 dimensional array from the start to the end. The array is 64KB (of 8 byte floats) and we assume the cache is 16KB:

[1,2,3.....2047,2048,2049,......4095,4096,4097.......8191]
 ^     ^        ^                    ^      ^
Cache  Cache    Cache                Cache  Cache
miss   hit      miss                 miss   hit

We are assuming that the CPU is implementing a simple caching system without any heuristics. When your code accesses the first element there will likely be a cache miss. However for the next 2048 elements the data will likely be in the cache already and calculations will progress faster. But for element 2048 we run out of data that is in the cache so that element will be slower to acccess while the value is fetched from main memory and the cache updated. Element 2049 will be relatively quick again until we hit element 4096. So there will be 4 cache misses as you process the array in sequential order. What happens when you process the array out of order? If you processed element 0, first then element 2048, then back to element 1? Well that could end up having a cache miss for every access - 2048 times slower!

So it is best to access elements that are next to each other.

The 2 dimensional case

Things get a bit more complex with the multidimensional case. By default numpy arrays are laid out in the “C” order. C uses row-major order, for an array with shape=(4, 5) this means:

array[0, 0], array[0, 1], array[0, 2], array[0, 3], array[1, 0], array[1, 1], array[1, 2].....

As you can see the last axis is grouped together in the array. So it makes the most sense to loop through each row so you are always accessing the next element:

for i in range(data.shape[0]):
    for j in range(data.shape[1]):
        sum += data[i, j]

If you were to reverse the order of these loops you would be ‘jumping’ about in memory by data.shape[1] - not a huge deal with small arrays but with larger ones this slowdown could be significant.

The multidimensional case and transposing

This only gets more important if you have arrays with many dimensions. Always loop through your array so the tightest loop is through the last axis, the second tightest loop through the second to last axis etc (assuming C layout).

This is how numpy does it internally.

But what happens when your algorithm requires visiting the axes in a different order? This is where numpy.transpose comes in. It is often worthwhile transposing the order of your axes first so that they are in the most optimal order for your algorithm. But won’t this create extra overhead? It does, but it is not as bad as you might think. Firstly the input array is read by numpy in a sequential order. Secondly, writing back to memory appears to happen “in the background” so usually by the time you are reading the values they have been written to memory and the CPU does not have to wait. Often there are fewer cache misses in total doing this. Will this help your particular use case? You’d have to benchmark the transpose and not transposed case and see what is faster. For big arrays it usually is faster.

Conclusion

Not all locations in memory take equal time to access. Because of the cache, most times it is the very next element from the previous one that will be fastest. You will need to think about this especially for multidimensional arrays. Note that processing using threads reduces predictability and often takes more CPU cycles for the same result. This is why Numba’s prange usually does not give expected speedups. Transposing your input array so it is in the correct order is often leads to a significant speedup, but the only way to be sure is to run timings for the various approaches.

Dealing with numpy arrays that have missing data is a challenge. When calculating statistics on them or otherwise processing data in them you need to skip elements that are missing - often they are some very high (or low) nodata value and will skew the result.

Using NaNs

In floating point data, there is a special number: NaN (“not a number”) that you can use to signify that this value can’t be processed for whatever reason. Numpy has a collection of functions (they all start with “nan”) that ignore any NaNs in your data. You an also test for individual elements being NaN with numpy.isnan. However, what happens if you are dealing with integer data? Well, setting integer elements to NaN fails. The alternative is to convert the whole integer array you are dealing with to float and back again. This adds time and memory use. Also performing operations on floats is slower than on integers.

Numpy Masked Arrays

There is an alternative - numpy masked arrays. These are arrays that also store a mask of where data is not valid:

>>> import numpy
>>> x = numpy.array([[1, 5, 9999], [-980, 7, 9], [11, 61, 9923]])
>>> mx = numpy.ma.masked_array(x, mask=[[False, False, True], [True, False, False], [False, False, True]])
>>> mx
masked_array(
  data=[[1, 5, --],
        [--, 7, 9],
        [11, 61, --]],
  mask=[[False, False,  True],
        [ True, False, False],
        [False, False,  True]],
  fill_value=999999)

It is important to note the mask parameter is True where the data is masked. Masked arrays support many of the numpy methods:

>>> mx.max()
np.int64(61)

Note how the masked out values aren’t used in the calculations. You can also turn a masked array back into a normal array using the filled method:

>>> mx.filled(-99)
array([[  1,   5, -99],
       [-99,   7,   9],
       [ 11,  61, -99]])

If your data has a single value that represents “no data” then you can pass this in to the masked_values function to create a masked array where that value is masked:

x = numpy.array([1, -99, 23, 78, -99])
>>> mx = numpy.ma.masked_values(x, -99)
>>> mx
masked_array(data=[1, --, 23, 78, --],
             mask=[False,  True, False, False,  True],
       fill_value=-99)

Notes on using Numba with masked arrays

Numba doesn’t know anything about masked arrays - you get an Unsupported array type: numpy.ma.MaskedArray error when you pass one in to a Numba function. However, a masked array is made up of two normal arrays: the data and the mask. You can pass these in separately to a Numba function:

@njit
def docalc(data, mask):
    tot = 0
    for x in range(data.shape[0]):
        if not mask[x]:
            tot += data[x]
    return tot
    
x = numpy.array([1, -99, 23, 78, -99])
mx = numpy.ma.masked_values(x, -99)
result = docalc(mx.data, mx.mask)

However, if your data just has a single “no data” value it may be easier just to pass this value in and compare each element to it instead of using masked arrays.

Conclusion

Numba masked arrays can be a useful tool when dealing with missing data. They provide a lighter weight alternative to conversion to float arrays and setting NaN.

Several of the packages available from the UBARSC group are most easily installed using the conda command. Some users are more familiar with using the pip command to download & install Python packages, and have wondered why we do not make our packages available in this way. What follows is a discussion of the major issues and reasons for this choice. Note that this is not a tutorial on the usage of either system.

What Are Pip, Conda, PyPI and Conda-forge

The pip command is a tool for installing Python packages. It understands the conventions used for distributing Python packages, and can be used to perform the installation directly from a source repository or .tar.gz distribution file. By default it will install things from the Python Package Index website (PyPI), handling all downloading directly. It understands when other Python packages are needed as dependencies for the target package, and will download and install those as well.

The conda command is a package manager with a much broader scope. It can install anything which has been packaged up for conda-based installation, and has no particular bias towards Python packages or anything else.

The PyPI website is a large repository for distributing Python packages, aimed at making it easy to download and install anywhere. It distributes both pure Python source, and pre-built “wheel” files (.whl) which might include compiled C-extension modules. It is the default source for packages to be installed by the pip command.

The conda ecosystem was developed by a private company, Anaconda (formerly Continuum Analytics), which sells (among other things) their services for software packaging and distribution. They have worked collaboratively with the open source community to maintain the conda tool and specifications, and host alternative community distribution channels. The major such alternative channel is conda-forge, and it is through this channel that we at UBARSC distribute our software. The conda-forge channel includes distributions for a wide range of data science related software, including binary distributions for many tools and libraries unrelated to Python.

Pip/PyPI vs conda/conda-forge

The PyPI repository is fine for distributing packages which are pure Python. All that is required is the Python source code, and some small text configuration files. However, distributing packages which are written in something other than Python is not really supported (although it can be kludged around, in some cases). A good example would be the KEA file format. The KEA library is written in C++, and while there is also a Python binding, this is largely separate from the library code itself. PyPI and the pip command don’t really have any idea what to do with all that.

Another more important example is the GDAL library. This is also written in C++, with an optional Python binding, and again, PyPI/pip has little to offer.

For example, RIOS depends very heavily on GDAL. So, if we were to distribute RIOS though PyPI, we ought to include this dependency in its requirements.txt file. However, if we did that, then pip would try to install GDAL from PyPI. The GDAL Python bindings are present in PyPI, but this is just the bindings, not the library itself, and so GDAL would still not be installed. If we do not specify the dependency, then the situation would be just as bad, and with little guidance for the novice user.

So, one would have to rely on something like conda to install GDAL itself.

The GDAL bindings also require numpy. However, numpy is also fully available from PyPI. This means that, depending on what else is installed, it is possible to have a numpy from PyPI and a GDAL from conda, and they may well be compiled in ways which are incompatible at the binary level, resulting in serious conflicts at run time.

Things can become even more complicated when one or more packages installed from PyPI include their own copies of binaries from other libraries. A good example of this is rasterio, which bundles its own copy of the GDAL binaries. These may be compiled with different options than any already installed version of GDAL, potentially causing all kinds of headaches.

To avoid much of this complexity, we have concluded that it is simpler just to rely on conda to distribute our packages, and everything on which they depend.

There is a great deal of discussion on some of these points spread all over the Internet. This article from Anaconda provides a brief overview (although slightly out of date).

As discussed in that article, conda and pip can work moderately well together, and this is useful for adding in small pure Python packages which are not available in conda. However, one does need to watch for pip’s tendency to try to install its own dependencies, and upgrade things which conda has already installed, resulting in incompatible combinations of binaries. For this reason it is strongly recommended that any package (including dependencies) which contains compiled binary files should come from conda-forge, and only pure Python should be installed from PyPI. Please exercise strong caution when combining them, including

1. being aware of the dependencies of the package you are about to install
2. watching the install as it runs, to see whether it tells you what else it is installing

This applies even when the target package being installed is from local source rather than PyPI. The ideal output of a pip installation should say that it installed the package itself, but that all dependencies were already satisfied. If pip starts over-writing things, you will probably need to discard that conda environment and start again.

Conclusion

In short, we recommend that you build working environments using the packages on conda-forge, including our own packages such as RIOS, PyShepSeg, TuiView, etc., and only use pip/PyPI to install small, pure Python packages, and only if they are not available via conda-forge. Of course, more sophisticated users have other options, such as building things directly from source, but we assume that such users are sophisticated enough to know what they are doing.

RIOS and PyShepSeg have support for Concurrency using AWS. RIOS supports ECS and Fargate clusters and PyShepSeg currently only supports Fargate clusters. These packages attempt to clean up any resources they create. However there may be situations, such as termination of the main script or software error, where these resources are not terminated. RIOS Reaper is a tool that monitors resources created by RIOS and PyShepSeg and notifies a user via email about anything that looks like it should be terminated. It does this by looking for resource tags that RIOS and PyShepSeg add as they are creating resources. It then checks EC2 instances for idle CPU and ECS and Fargate clusters that are stopped. Currently, RIOS Reaper does not terminate any EC2 instances or Fargate clusters. We feel it is safer for a user to delete these in case they are still in use.

Installation

RIOS Reaper is a Lambda that runs once a day. It is based on AWS SAM and this needs to be installed first. Any machine that is connected to the internet and logged into AWS should be suitable for running the deployment, however we have only tested on Linux. Further instructions are in the README. Note will need to be logged onto your AWS account on the command line with sufficient parameters to install a Lambda.

Since this job runs on a Lambda at a nominated time there is no extra machine to be provisioned and you’ll only be paying for when the Lambda is running. AWS takes care of running Lambda code, you just need to provide the function.

Information about your account and VPC is required. This is passed in via environment variables:

• AWS_PROFILE - the name of the profile you are running under, or default
• VPC_ID - the id of the VPC you want the Lambda to run under
• SUBNET_IDS - a comma separated list of subnet ids the Lambda is to run within
• EMAIL - the email address you want notifications to be sent to

You may also want to make changes to template.yaml which is the CloudFormation template used by AWS SAM. Check that the Architecture parameter is correct for your configuration (also in samconfig.toml) and that the Schedule setting is correct for when you want the check to be run. The region can be set in samconfig.toml.

Once this information is set, you can test running the Lambda locally:

./test-deploy.py

Once you are happy, deploy the Lambda like this:

./test-deploy -m deployed

To remove the Lambda:

sam delete

Conclusion

RIOS Reaper is a relatively easy way to check that there are no “zombie” jobs running that are costing you money. You will get an email once a day containing details of instances and clusters that may need to be checked and terminated.

Although TuiView is a useful program on its own, sometimes it may be useful to embed some of its functionality within another script. In this way you can build a customised image viewer. Since TuiView is just a Python module you can access it from any Python script. However, a certain amount of Qt/PySide knowledge is required for building user interfaces. This tutorial is based on the TuiView wiki. Documentation for the TuiView internals can be found at the TuiView Developer Documentation.

Creating a TuiView widget within your own window

This is the simplest scenario. You have your own window, but you want one of the widgets to be a TuiView map:

Below is the source for this application:

import sys
from PySide6.QtWidgets import QApplication, QWidget, QVBoxLayout
from tuiview import viewerwidget, viewerstretch
from osgeo import gdal

gdal.UseExceptions()

app = QApplication(sys.argv)

ds = gdal.Open('a.tif')
stretch = viewerstretch.ViewerStretch()
stretch.setRGB()
stretch.setBands([4, 3, 2])
stretch.setStdDevStretch()

w = QWidget()

tmap = viewerwidget.ViewerWidget(w)
w.show()

# Note that TuiView expects the widget to be shown before
# adding a layer
tmap.addRasterLayer(ds, stretch)

layout = QVBoxLayout()
layout.addWidget(tmap)
w.setLayout(layout)
w.resize(250, 150)
w.move(300, 300)
w.setWindowTitle('Simple')

app.exec()

Creating new viewers

If you are happy with the way the TuiView windows look, but you just want to drive them programmatically, the recommended way to do this is to use the GeolinkedViewers class:

Below is the source for this application:

import sys

from PySide6.QtWidgets import QApplication

from tuiview import geolinkedviewers
from tuiview.viewerwidget import GeolinkInfo

app = QApplication(sys.argv)

glviewers = geolinkedviewers.GeolinkedViewers()
viewer1 = glviewers.newViewer('a.tif')
viewer2 = glviewers.newViewer('b.kea')

# The first parameter is the 'id' of the sender; 
# set to 0 if not sent from inside TuiView.
# Then pass Easting, Northing and meters per pixel as zoom factor
obj = GeolinkInfo(0, 1976486, -3144006, 100)
glviewers.onMove(obj)

app.exec()

Conclusion

Most aspects of TuiView can be automated from a Python script. This means you can create your own custom applications that re-use TuiView functionality without having to reinvent the wheel.

This script reads one or more Python files and writes corresponding Markdown files containing their docstrings, formatted in a neat and readable form.

It is intended as a simple alternative to using tools like Sphinx, ReadTheDocs, or GitHub Pages, when the complexity and overhead of those tools is not justified.

The advantage of writing Markdown is that it can be viewed directly from a Github repository. Thus, Github can easily serve the documentation pages for a project, without setting up any extra systems.

Although the UBARSC projects rely heavily on the Raster Attribute Table (RAT) functionality within GDAL, we accept that for some situations this can be problematic:

1. RATs with many columns. Eventually this becomes unwieldy. It would be nice if columns could be grouped in some way and just access what would need.
2. Updating RATs is problematic for files on S3. They must be copied locally, updated then copied back. It would be nice if you could just write a new column to the file on S3.

We have identified the Zarr Format as a useful alternative and are working to support RATs in Zarr files alongside RATs in normal GDAL files. Among the useful features of .zarr files, they can be updated on S3 directly.

What is changing?

If you are currently using RATs in GDAL files (like KEA and HFA) and you are happy, then nothing will change.

Introducing RatZarr

RatZarr is a new UBARSC project that creates a RAT-like interface to a Zarr file. Note that it doesn’t work with any arbitrary .zarr file - only ones created by RatZarr itself. Most users won’t use this project directly - only via RIOS and/or pyshepseg. Note that we aren’t talking about saving imagery into a .zarr file - just the RAT.

RIOS 2.0.9

RIOS 2.0.9 has been released with support for .zarr files (via RatZarr, if installed) in ratapplier. The full list of changes can be found here. An example of reading and writing data to/from a .zarr file (alongside a KEA file) is below:

from rios import ratapplier
inRats = ratapplier.RatAssociations()
outRats = ratapplier.RatAssociations()

inRats.vegclass = ratapplier.RatHandle('vegclass.kea')
inRats.heightclass = ratapplier.RatZarrHandle('heightclass.zarr')
outRats.vegclass = ratapplier.RatHandle('vegclass.kea')
outRats.heightclass = ratapplier.RatZarrHandle('heightclass.zarr')

ratapplier.apply(myFunc, inRats, outRats)

def myFunc(info, inputs, outputs):
    outputs.vegclass.colSum = inputs.vegclass.col1 + inputs.vegclass.col2
    outputs.heightclass.colSum = inputs.heightclass.col1 + inputs.heightclass.col2

Note that .zarr files can be accessed on S3 using the s3://bucket/path/to/file.zarr syntax.

PyShepSeg 2.0.5

PyShepSeg 2.0.5 has been released, also with support for RATs in .zarr files (via RatZarr). The full list of changes can be found here. You would use this functionality if you wanted to save statistics to a .zarr file, like this:

from pyshepseg import tiling
segResult = tiling.calcPerSegmentStatsTiled('veginfo.kea', 1, 
    'segment.kea', statsSelection=[('Band1_Mean', 'Mean')],
    outFile='s3://bucket/veg.zarr', outFileIsZarr=True)

You then would presumably do more processing with the data in the .zarr file with rios.ratapplier as detailed above.

Conclusion

Reading and writing RAT columns to and from .zarr files is a very useful feature, especially for those working in cloud compute environment. New releases of RIOS and PyShepSeg allow users to save their data in .zarr files and any feedback is appreciated. Future development will include support in TuiView.