ILWISPy

Notebook prepared by Ben Maathuis, Martin Schouwenburg and Bas Retsios. ITC-University of Twente, Enschede. The Netherlands

ILWISPy is a Python library for processing, accessing and writing geodata. The library handles raster data (2D/3D), vector data and tabular data. The core/backend of the new ILWISPy is a general purpose GIS library written in C++. Its modular architecture (supported by Qt) allows programmers to write native plugins (C++) to extend the libraries interoperability or processing functionality. The user can import the whole ILWIS functionality into Python, which facilitates writing high level scripts for multi source geo processing in a standard scripting environment. ILWISPy is running on different platforms, e.g. on Linux as well as on Windows based operating systems, for Python versions 3.6 and higher. ILWISPy can be directly imported in Python, but can also be used within a Jupyter Notebook, offering the capability of markdown text with explanations and subsequent coding fields to execute certain operations. This tutorial therefore provides a brief introduction into the ILWISPy API as a plugin in Python, executed within a Jupyter Notebook.

This notebook will highlight how to obtain information on the various operations available within the library and how to use the fuctionality provided using the ILWISPy library. For visualization of the results you can use ILWIS-386 desktop version, available at: https://filetransfer.itc.nl/pub/52n/ILWIS386/Software/. Special attention in this notebook is given to the interoperability of ILWISPy with a number of common used Python libraries.

All ILWISPy releases are adhering to the GNU GENERAL PUBLIC LICENSE - GPL3 (https://www.gnu.org/licenses/gpl-3.0.en.html). The ILWISPy development team is consisting of Bas Retsios, Martin Schouwenburg, Willem Nieuwenhuis, Lichun Wang and Ben Maathuis, from the Faculty ITC, University of Twente, Enschede, The Netherlands. For questions or further information don't hesitate to contact the ILWISPy team members Bas Restsios or Ben Maathuis at v.retsios@utwente.nl / b.h.p.maathuis@utwente.nl respectively.

Why ILWISPy

The way currently remote sensing images, derived products and data is provided has changed the ways how this information is being processed. Traditionally a stand-alone system was used with a number of software tools to process / analyse the data. Nowadays information from the cloud is accessed and processing is done using tools provided by these cloud based service providers. Once analysis has been conducted, the final products are downloaded to integrate these with existing, local data sources for final analysis. In short some of the advantages bulletwise:

Within this notebook an introduction is provided to a number of topics:

Additional installation instruction and updated information can be obtained at: https://www.itc.nl/about-itc/organization/scientific-departments/water-resources/software-tools-models/ilwis3-and-toolbox-plugins/ilwispy-getting-started/. Special attention has been given to the use of Python to show the seamless synergy with ILWISPy functionality to process geo spatial information. For visualization of the data used in the notebook or of the processing results obtained use can be made of ILWIS 3.8.6.

ILWISPy under Python - Jupyter Notebook

Installation instructions when using Windows OS:

First steps using ILWIS under Python requires installation of ILWISPy as a site-package. If working using a Linux Ubuntu Operating System, a Linux - ILWISPy version for Python 3.8 is available. For Windows Operating Systems various Python based versions are supported, e.g. Python 3.6 up to Python 3.11.

The wheels can be downloaded from 'https://filetransfer.itc.nl/pub/52n/ilwis_py/wheels/ ', note the "whatsnew.txt" file, check the information provided, subsequently select the folder related to your operating system (windows or linux) and navigate to the most recent dated folder to download the required wheel for your Python version.

Install the wheel corresponding to your Python version, as example below the ILWISPy version for Python39 (64bits) is used. Using the command prompt from the Python39 folder, type the following command:

"python -m pip install Ilwispy* - *.Whl"

note: - is representing the wheel version used, here for python39, 64 bits, e.g. ilwis-1.0.20230330-cp39-cp39-win_amd64.whl. If you have installed a previous version of ILWISPy, you can use the procedure described above to reinstall a more recent version, the older version will be removed and the new version will be installed.

It is assumed furthermore that the Jupyter Notebook site package has already been installed.

To check if your ILWISPy installation was successful, use the Command Shell, navigate to your Python folder, start Python and within the Python Command Shell execute the following expressions:

Use quit() or Ctrl-Z plus Return to exit.

Now you are ready to explore ILWISPy.

Quick overview

The figure / table below provides an overview of geospatial data and the relationship between them as used in ILWISPy. For example raster data, here called a 'RasterCoverage' requires next to the actual raster data, a GeoReference and a DataDefinition. Same is applicable to vector data, although other elements might be required. Review the figure below for the linkage of the various geospatial elements, their description is given in the table.

IP_spatial.png

Geo Spatial Element Description
RasterCoverage Gridded raster data in 2D or 3D. Numeric values or thematic data
FeatureCoverage Vector data in points, lines and polygons with attribute data
Table Tabular data for numeric, string and thematic data
DataDefinition Defines what type of data is used and what valid entries (range) are available
Domain The datatype of the data. Can numeric (various), string, identifier/thematic (various), date/time etc...
Range Describes the valid range of values that are allowed for certain DataDefinition. e.g the number 0 to 100. For numeric data this includes the resolution of the numbers.
GeoReference The Geometry of a RasterCoverage. What is the connection of the pixels in the grid and the real world coordinates.
Size Simple object to hold X/Y/Z size.
Envelope Object to define the area (rectangular) in real world coordinates of an object
CoordinateSystem Definition of a real world reference system to define locations.
ColumnDefinition Simple combination of a datadefinition and a columnname. Is the same for features and tables though in the former case we often use the name attribute.

At your disposal are a set of functions, which are called 'Operations' and provide access to all the processing capability available within the ILWISPy library able to handle and manipulate the various geo spatial elements indicated above.

Getting familiar with the ILWISPy Library

First step using ILWISPy under Python is loading the required libraries, including some of the required python libraries, like os, for the use of operating system dependent functionality, numerical python and matplotlib, a comprehensive library for creating static, animated, and interactive visualizations. Next to ILWISPy use is made of this added functionality.

To check the ILWIS version activated as Python site package

Different types of operations

Within the ILWISPy Library a collection of operations are at the users disposal. There are operations that provide information about the object and there are operations that perform a manipulation on source data provided and create new output.

A number of operations are a function of the object. Examples related to raster objects are: array2raster, pix2value, statistics and size, etc. Some of them will be further introduced below. These operations can be used to retrieve metadata, these 'object.operations' never modify the object itself.

Working with geo-spatial data using ILWISPy - data mining and metadata handling

Let us see if we can access some raster data and get familiar with some of the objects mentioned above. If we work with disk based data it is often practical to set your working catalog. The use of path names to find your data is no longer required.

Get a listing of the catalog - folder for a specific item, below a listing is provided for all raster layers, (.it.RASTER) can be replaced by ilwis type such as POINT, FEATURE, POLYGON, etc.

This shows all rasters in the current working catalog. ilwis.catalogItems() will show all objects in their roles. This sometimes can be a bit confusing sometimes as an object can appear multiple times in the list in different roles. E.g. a GTiff file can be a raster, obviously, but it also can be (used as) a Georeference as it contains all information to function as a GeoReference. See what happens

So it it often wiser to use a filter to know exactly what you are viewing/searching. Let's now turn to how to access the data. We start with data in a raster or regular grid format.

Raster data

It is advised to familiarize yourself, before you contine, with the sample data used in this notebook. You can use ILWIS 3.8.6 (available at: https://filetransfer.itc.nl/pub/52n/ILWIS386/ - in the folder 'Software'). If not familiar with with ILWIS check also the folder 'Tutorial' for additional information on installation and using the desktop version.

Working with single layer raster data

Having set the location where to find the data (workingCatalog) we don't need a full path any longer and can open the raster coverage directly, even if the raster file is not in a native ILWIS format. Here the initial file format is GeoTif.The raster data is containing NDVI values in a range from 1 to 255 - represented from light yellow to dark green respectively, the water and areas outside Africa are assigned as no-data and are assigned 'white' in the figure below. Further below you will also create raster visualization, here a picture of the raster is given for your reference.

NDVI image over Africa NDVI.PNG

Loading the raster data

Check the raster dimensions

Apparently it is a raster of X size 1152 pixels, Y size 1152 pixels and contains only one band, in this case the NDVI values of Africa and part of the Middle East. An ILWISPy raster may contain many bands and is still considered one rastercoverage.

This is the size of the raster in real world coordinates. Order of corner coordinates is Min X, Min Y, Max X, Max Y. Obviously not latlon so lets try to find out what kind of coordinate system is used.

Apparently it uses a projection that is 'Albers Conic Equal Area' and 'Clarke 1866' as ellipsoid. Another well known definition is if we try to find its proj4 definition. PROJ is a library for performing conversions between cartographic projections.

Though there is much more to say about the exact geometry of the data lets now turn to the data itself. First let us examine what type of data of the raster. As was shown in the diagram above the type of the data is encapsulated by the datadefinition. So let's start there

Every datadefinition contains an attribute that is called 'Domain'. It encapsulates the type of data that is describes by this definition. Here the data type is 'image' or a Byte / Unsigned Integer 8 data type . Another common data type / domain is 'value'. ILWISPy internally uses a float-64 data type in memory, but upon data storage on disk the data type and precision can be defined by the user.

Let's have a look at the data contained in the raster into more detail, e.g. by creating the raster statisics. First create a histogram listing. Note that the binning interval used is 2

Note the last tuple (-1e+308, 803014.0) which actually indicates if 'no data' is appearing in the image, which is clearly the case here as there are quite a lot of pixels which are assigned no-data (803014.0). Below we create additional statistical objects that contains some of the other basic descriptive statistical information for this raster.

To plot the histogram, the commonly used Python library 'Matplotlib' can be used. The library was already loaded and a plot can be created. Note the last tuple (containing the no-data) is not taken into consideration when creating the plot

It is also possible to retrieve the value at a certain X, Y location, taking into consideration that the row/column index will start from 0,0. Note that when reviewing your result using the map in a desktop GIS the row / column numbers likely will start at 1,1 for the first pixel in the first row. Therefore to check manually the value in the map you will have to correct for this offset, so the actual column / row number is 501/701!

This gives the value of pixel at location (x,y) 500,700. Note that pixel locations start at 0,0. The function is fairly expensive (in terms of performance) but useful when inspecting individual locations in a raster. A better way if one wants to do some (simple) processing is the 'PixelIterator'.

Loading an ILWIS raster file (extension *.mpr) containing elevation information and retrieving the meta data and descriptive statistics

It is also possible to create a definition in order to retrieve and display the main descriptive statistics and plot the histogram as this type of information is required more often, only for different rasters

To execute the definition the required raster input layer has to be selected

Plot the raster data using the Python library

ILWISPy is fully compatible with Python, this also is applicable to raster data. Below the DEM in ILWIS (.mpr) format is transformed into a Numpy array and is visualized using Imshow(), which allows us to plot an image from a 2-dimensional numpy array

Using multi layer raster data

ILWISPy makes no difference between single band files and multiband files. In fact the single band file is just a special case of multiband, having one band only. Below a multiband raster in ILWIS format is loaded, note the extensions of the file below, '.mpl' which is refering to a 'map list'.

A multi raster containing 12 bands has been loaded. Everything works the same apart from the fact that we now have a Z component in the pixels, the band index starts at '0'. A value in the fourth band (remember the counting starts at 0, so 3 is the fourth band) can be retrieved in a simalar way as a single band raster, now including the layer number.

Retrieve a certain value contained in any of the raster layers. The iterator now works over all the bands (processing might take a bit longer).

Apparently in band 10 and 11 (representing layer 11 and 12 in the multiband raster!) the value 2.40 can be found. Later when discussing 'operations' you will see that processing single or multiband rasters is identical.

Another example to obtain information of the same pixel but for all layers in the multibandRaster loaded is given below. Check the values in the list obtained and compare with the result above.

Create a simple plot showing the Monthly Temperature over the pixel location selected

Creating new ILWIS raster data and the exchange between ILWISPy and Python - Numpy

Sometimes you need to create an empty raster in order to fill it with data comming from another source (e.g. from a Numpy array). So how do we do that? Looking at the diagram in the beginning we need a GeoReference (geometry), a DataDefinition (Data type) and the actual raster data.

The first line creates a GeoReference with a certain coordinate system (epsg:4326 = LatLon), certain bounds and size in pixels. The second line creates the DataDefinition. It has a numeric value domain ('code=value') and may contain integer numbers between 0 to 1000. Attach these elements to an emptry RasterCoverage and you have obtained a functional raster.

Next we are creating some synthetic raster data. This is just an example so its rather simple, adding the value of 1 to each next raster pixel on a line by line basis, starting from the top left (first pixel) till the bottom right (last pixel) for the array size defined, using a float64 data type

So how do we get data into this RasterCoverage? There are several methods to do so, but the 'array2raster' is used here to copy the data to the newly created raster 'rcNew'

Get dimensions of this newly created ILWIS array

We can also convert the ILWIS raster array back to a numpy array, so it can be visualized using the matplotlib function imshow(). First the ILWIS raster is transformed into a numpy array

Note the order of Z,Y, X. Note that here Z is only 1 layer, same could be applicable to a multi image stack

To display the numpy array use can be made of the matplotlib function imshow()

Check value of last pixel just to show that data has been properly processed, do you agree with the value retrieved?

Writing a variable to a file

Though there are of course many details and options that are left unexplained for the moment there is one point that needs to be addressed. How can I store data in an appropriate data type, having the required precision. Below the 'ilwis.do' operation 'setvaluerange' has been included, prior to saving the data. There will be additional information on the 'ilwis.do' operations available in the ILWISPy Library further below, here it is sufficient to review the parameters that are required, like the minumum, maximum values of the data range as well as precision for the data to be stored.

Here we store the just created data in the 'ilwis3' 'map' format ( rasters) with the name 'mydata'. If you don't include a path it will go to the working catalog else it will go to the path specified. The format name is a name that is relevant in the context of the 'provider'. The provider is the last parameter. A provider usually represents an underlying library that does the actual storing to the requested format.

If you want to visualize your raster data in the ILWIS desktop package, e.g. ILWIS 3.8.6, the statement above can be further simplified

The statement below does the same thing but this time it uses the gdal library to store the output raster in a GeoTiff format (for GDAL (short)format names see: https://gdal.org/drivers/raster/index.html and https://gdal.org/drivers/vector/index.html).

Vector Data

Now have a look at how ILWISPy deals with vector data. Vector data which is used in features is loaded in the same way as raster data. It can come from disk based files or from a Postgress database. In this case we use a shape file available in the tutorial data folder.

This loads a feature coverage containing towns (as points) in Ethiopia. So lets examine the data

Apparently we have 934 points in the coverage. Note that FeatureCoverages are not restricted to one vector type. One could have points and polygons in the same coverage. Though in this case we only have points. We could query.

As with raster the geometry of the coverage in world coordinates is determined by its coordinate system

Features can have multiple data attributes

Lots of output. But basically this tells you what attributes are available and what their definition is. As with rasters there is also a statistics object (for numerical attributes).

So what about the features themselves? How dow we reach them? Like with RasterCoverages we use an iterator for this. With the iterator we can access individual features. On top of that there is a '[]' operator to access the value(s) of the individual attributes of the features.

An alternative way of accessing the data is viewing the coverage as a indexable data array and using this model to access its data.

Because we can iterate over the features in a feature coverage the normal loops in Python will work too

Changing attribute values is also possible by using the setAttribute method of a feature.

Note that the attributes can also be retrieved by something that is called the attribute table. Basically it is a representation of all the attributes of the features as a table with records for features and columns for attributes.

In most respects it's a normal table but the backend is the featurecoverage and not a copy of the data of the featurecoverage. So changes on the table, changes the coverage.

The one attribute we didn't talk about is the geometry of a feature. How is defined what the positional structure of a feature is? The shape of a feature is defined by its geometry.

The WKT method gives the (OGC) WKT definition of the shape/geometry of the feature in terms of the coordinate system attached to the coverage. Obviously for lines and polygons this is a much(much) longer string. One can also set a geometry in this way. (ignore the rounding stuff far beyond the important part of the numbers)

How about creating a totally empty FeatureCoverage and adding features to it. In the example below we create a new FeatureCoverage from the old one with only those cities over 100000 population.

In the newFeatureFrom method, the second parameter (coordinate system), facilitates automatic coordinate transformation if the source map (fc) doesn't match the coordinate system of the target map (fcNew). In this case they match.

We can also create and fill a completely new FeatureCoverage.

We now have an empty coverage with two attributes : 'Name' and 'Price'. Next we add some values to it and check if the results are what we expect.

As before with the RasterCoverage we skip some topics for brevity. The last point I want to touch is storage of your FeatureCoverage. The method is identical to that for rasters

Working with a Feature Coverage - Polygon map

Retrieve the classes of the polygon map

Fetching the attribute details for the first feature from the polygon map

Retrieving features and storing the selected features as a new shape file

Tables

The last of the primary data objects is the table. As with the other data objects it's source may come from file based data (ilwis3 tbt, excel sheets) or from a postgres database.

This is similar to how FeatureCoverage retrieved its attribute structure. For the rest there is not much to say about the meta-data of a table. The data itself can be access/changed via the 'cell' or 'record' methods.

The table allows some simple selections via the select method. It returns the record numbers for which the conditions holds true. These numbers can then be used the to access the data itself.

Changing the data happens through the setCell methods. The first parameter is always the column (either by index or by name, then the record number and then the new value.

Operations - ilwis.do

ILWISPy allows using simple (single line) syntax, which otherwise under Python would require many lines of code, to conduct common GIS and RS operations. Their primary purpose is to be used in the processing operations that ILWISPy supports. The exact number of operations can vary, e.g. depending on which plugin versions are used. A listing of the current available processing operations can be obtained by:

Note that some operations appear double(or more) because operations with same name exist but with different numbers of parameters. Before delving deeper into how we work with operations lets first explore how we can get some extra meta-data from the operations to understand how they work.

For a more appropriate alphabetical listing on the available operations, including information on the parameters required some additional Python code can be applied

To get additional information on an operation, the input and output parameters and well as a description of the operation the following code can be applied. Here as example the 'resample' operation is used. To extract meta-data there is the function 'operationMetaData'. It can extract meta-data per operations. Its base call extract the syntax of the operation ( later we will see how to use that).

At least now you have a rough idea of what this operation needs. A raster, a target georef and an interpolation method. There are some other properties we can query. We can also query the individual parameters to see what they mean.

Hopefully this metadata gives enough of a clue how the operation works if you are not familiar with the operation. Now we are going to tackle the actual execution of the operation. All operations use the 'do' method. Lets try the resample method. The 'ndvi.tif' map is a continental map in the Albers projection and could easily be transformed into the 'LatLon'coordinate system. Note we use a previously unmentioned method of CoordinateSystem to convert an envelope to an envelope in another CoordinateSystem : 'convertEnvelope'

Now we can do the actual operation and fill in the parameters according to the syntax

To get additional information on an operation, the input and output parameters and well as a description of the operation the following code can be applied. Here as example the 'timesat' - time series filtering operation is used. Basically this all there is to be said about operations. Every operation works the same, no exceptions :

outputdata = ilwis.do(operationname, inputdata, parameters)

Next we examine a few operations that might be of importance though they all use the same method('do') to work. The main difference is that one of their parameters is like a mini script, allowing a very flexible kind of operation. Lets first have a look at the 'mapcalc' operation. This operation can do mathematical calculations on numeric rasters. Note that the rasters need to have the same Geometrty (GeoReference). E.g.

In this simple statement every value in the raster gets 100 added and the result is written to the new raster. The @1 points to the first input parameter that follows the expression. @2, @3 etc... would point to the subsequent parameters (max 6)). The expressions might get more complex with inbuild operations, brackets and multiple input maps.

An import role in the 'mapcalc' semantics is played by the 'iff' in build function which can fork decisions based on the values in the raster

In this case all values above 40 become undefined and the rest becomes 1. mapcalc, nor any other operation, cares about how many bands a raster has. The operation applies to all bands

Obviously it takes longer but all the pixels for all bands in the multiband raster have been recaculated. The following built in operations are available : iff,sin,cos,tan,asin,acos,atan,log10,ln,exp,abs, ceil,floor,sq,sqrt,max,min,pow,not,xor,ifundef, ifnotundef

Next we look at the 'selection' operation for rasters. This operation 'selects' pixels in an input raster and write the selected pixels to the output raster. E.g

This cuts(selects) a section of the raster and puts it in a new raster

The same but this time only the pixels with value greater than 40 will be selected. Note that 'pixelvalue' is a special attribute that indicates the actual value of the pixel. If one uses a raster with an attribute map one also may use the attribute values.

Again the same but we define the selection as a range ( between 40 and 150). Note that the spatial contraint can also be an envelope ( so coordinate box) or a wkt definition of a polygon.

As usual there is a lot more to say but for an introduction to operations I think this suffices.

Geometric Operations

The geometries support a number of the operations as defined by OGC. e.g. The distance between two features

The operation supported are distance, getArea, getLength, buffer, convexHull, intersection, difference, symetricDifference and join. E.g.

The boolean intersection operation supported are within, contains, disjoint, touches, intersects, crosses, overlaps, equals, equalsExact, covers, coveredBy and relate as defined by the DE-9IM intersection matrix.