ILWIS Python API Tutorial

About

The core/backend of the new ILWIS Objects 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. Next to native plugins, ILWIS Objects also supports the well documented ILWIS-script, which allows ILWIS users to write processing scripts, instead of hardcoded plugin functionality. Finally the user can import the whole ILWIS functionality into Python, which allows him to write high level scripts for multi source geo processing in a standard scripting environment, possibly running on different platforms. This tutorial provides a brief introduction into the API ILWIS provides as a plugin in Python.

Download & Installation

The ILWIS Objects Python extension is only available for the 32 bits versions of Python 3.3, 3.4 and 3.5 on Microsoft Windows. We cannot support 64 bits Python 3.x yet, nor Linux or Mac, or Python 2.x.

1. Python 3.x (win32)

   1. Download the latest 32 bits Python 3 installer from https://www.python.org/downloads/ and start the installer wizard.
   2. IDLE (Python GUI) might be useful if you have no other Python IDE installed. IDLE is included in abovementioned installer, when choosing the default installation.

2. ILWIS Objects Python extension

   1. Download ILWISObjects-beta3.0-forPython3.x-32bit.exe that corresponds to your Python version from https://github.com/52North/IlwisCore/releases/tag/v1.0-beta3.0
   2. Follow the ILWIS Objects install wizard, you will be asked to select the Python 3.3, 3.4 or 3.5 installation directory of a 32 bits Python installation, to configure the extension. Note that this folder (usually C:\Python3x) is different than the folder where ILWIS Objects is installed (usually C:\Program Files (x86)\n52\ILWISObjects).

3. ILWISObjects-HelloWorld.zip
   Download from https://github.com/52North/IlwisCore/releases/tag/v1.0-beta3.0

Hello World

The ILWISObjects-HelloWorld.zip archive contains the Python script (helloworld.py), which we go through here, plus the ILWIS3 feature dataset of the world countries. Please unzip everything before going on.

Explaining the feature data set

The ILWIS3 data format includes multiple files containing either meta or binary data, just like the ESRI Shapefile format. The multi-polygon geometries representing the world's countries are stored in countries.mpa (meta data) and in countries.mpz# (binary data), whereas the feature attributes like name and population (pop_est) are stored in (attribute-)table files called countries.tbt (meta data) and countries.tb# (binary). In ILWIS each dataset is called by the meta data file name.

A first ILWIS-Objects-Hello-World-Python-script (helloworld.py):

import ilwis
ilwis.Engine.setWorkingCatalog("file:///D:/<working>/<directory>")
world = ilwis.FeatureCoverage("countries.mpa")
population_ranking = {}
for country in world:
  population_ranking[str(country["iso_a2"])]=int(country["pop_est"])
print(sorted(population_ranking.items(), key=lambda x: x[1]))

To run this script, please open it in IDLE (Python GUI), insert the path to the unzipped data folder into line 2 and hit F5 or <Run><Run module>.

This script should print the sorted list of tuples representing each country in a rank ordered by its population like:

[('CN', 1338612970), ('IN', 1166079220), ('US', 313973000), ...

Explaining the Code

In line 1 we import the ILWIS library (called ilwis).

In line 2 we set the current working catalog to a directory in a local file system, indicated by a URL containing of the protocol prefix file:// followed by a directory path /D:/<working>/<directory>/ to the place, where the countries.* files are. For a local file system on Windows, the directory can also be specified as D:\<working>\<directory>. As soon as the working catalog is set, ilwis scans it for relevant geo-data files, and probes for their metadata. The metadata is remembered in order to speed-up reading and processing later on. In order to recognize catalog items as usable geo-/meta- data, ilwis internally relies on so called connectors (plugins). Since the standard installation of ILWIS Objects comes with an ILWIS3- and a GDAL-connector, it will detect a large number of vector and raster formats, for example *.shp, *.kml and *.tif as well as their ILWIS3 equivalents *.mpa, *.mps, *.mpp and *.mpr.

Now we can start analysing the vector file countries.mpa as in line 3 we loaded it into a FeatureCoverage object, which now manages all the geo-information about the features contained in the file (in this case countries).

In line 5 our simple analysis uses the standard Python iterator interface (__iter__()) of FeatureCoverage to loop over each country in the world and stores (in line 6) the population attribute (country["pop_est"]) as a native Python int object into the dictionary population_ranking using the str typed ISO abbreviation attribute (country["iso_a2"])) as key.

Line 7 is a very pythonic way of sorting the dictionary population_ranking by converting it to a sorted list of tuples (.items()) using the second part of each tuple as the sort key (key=lambda x: x[1]) to finally print it to output.

IlwisObject

Most IlwisObjects allow to load data from files or other external sources but by default, ILWIS will also create a set of internal objects to be used on the fly. Also you can create a new IlwisObject like an empty FeatureCoverage just by calling its constructor and giving a new name:

fc = ilwis.FeatureCoverage()
fc.name("newFC")

You can later ask each IlwisObject about it’s name:

fc.name()
'newFC'

or change it to a different one:

fc.name("newName")

fc.name()
'newName'

This new FeatureCoverage can be edited by adding/changing features and attributes. Unless it is stored to disk, it will stay an internal (temporary) object until it is no longer required / runs out-of-scope. You can check if the object was stored:

fc.isInternal()
True

To store an IlwisObject to an external resource, use the .store(...) function:

fc.store("file:///D:/<dir>/countries.shp", "ESRI Shapefile", "gdal")

Depending on the connector's implementation the given URL can either point to a directory, a file or whatever writable resource. The lines above will cause the GDAL connector plugin to store the typical set of ESRI Shapefile format files (countries.shp, *.dbf, *.prj and *.shx) into D:/<dir>/. To store the same object to the ilwis3 format (ILWIS polygon .mpa files):

fc.store("file:///D:/<dir>/countries", "vectormap", "ilwis3")

Since version Beta 3.0 of the Python extension, it is possible to simplify the store():

fc.store("file:///D:/<dir>/countries.shp", "ESRI Shapefile") # will store using gdal
fc.store("file:///D:/<dir>/countries", "vectormap") # will store in ilwis3 format (because the configured ilwis3 connector recognizes "vectormap" as a valid format)
fc.store("file:///D:/<dir>/countries") # will use ilwis3 (our default) and automatically determine the format to store

Coverage

AttributeTable

FeatureCoverages and RasterCoverages allow access to their corresponding AttributeTables as separate IlwisObjects like this:

table = world.attributeTable()

But you can also create a stand-alone Table from a file, like this:

table = ilwis.Table("rainfall.tbt")

Using the Table class, attributes can be manipulated on a cell-by-cell-basis:

value = table.cell("iso_a2", 4)
table.setCell("iso_a2", 4, value.upper())

Here we transform the value of record number 4 in column "iso_a2" to an upper case string.

Furthermore it is possible to request the structure of the Table, like this:

recCount = table.recordCount()
colCount = table.columnCount()
columns = table.columns()  # ('column1','column2',...)

Where .recordCount() returns the number of records (rows) and .columnCount() the number of columns in the Table. In the last line, .columns() will return a tuple of column names (str). Instead of column names you can always also use its column index (unsigned integer).

!

Metadata that is not stored in the data header (depending on the format), might not be correct if the actual binary data was not yet (fully) loaded, so methods like recordCount(), featureCount() or envelope() may contain wrong data, until ILWIS is forced to load the binary data, by accessing the first record, feature or pixel.

column = table.column("iso_a2")
column = table.column(42)
columnIndex = table.columnIndex("iso_a2")

The first two statements return exactly the same data, which is the whole column as a Python Tuple. If you don’t know the index, you can call .columnIndex("iso_a2")to retrieve the desired index. Similarly you can access the data of a record by index like this:

record = table.record(4)

This way you will get a Tuple of the whole data of a record.

There are two ways of adding a new column to an attribute table. One way looks like this:

world.addColumn(coldef)

Here you just pass in a new ColumnDefinition. And the second way is:

world.addColumn("newCol", "value")

where the first argument defines the name and the second the DataType/Domain of the new column. Note that the addColumn function is applied to the FeatureCoverage object, thus not on the Table! This is to make explicit that the schema is to be altered at the FeatureCoverage. For stand-alone tables however (thus not attribute-tables), the addColumns is applied directly on the Table object:

table.addColumn(coldef)
table.addColumn("newCol", "value")

As mentioned, these addColumn() calls are not allowed on tables that are attribute-tables; for those, the table-schema can only be altered through the feature-coverage.

!	Trying to add columns to a table that is retrieved from a `FeatureCoverage` with `.attributeTable()` will provoke an error.

The codes of the following internal standard Domains can already be used for column definitions:

CODE	DESCRIPTION
text	strings
value	numeric values with variable step size fitting into a double (floating point)
count	positive integer values from zero to 999999999
integer	integer values
boolean	True or False
percentage	0.00% to 100.00% (step size 0.01)
distance	0.00 to 1000000000.00 (step size 0.01)
NDVI	-1.00 to 1.00 (step size 0.01)
time	ISO 8601 compatible times

CoordinateSystem

CoordinateSystem objects are either created from a Proj4 definition:

cs1 = ilwis.CoordinateSystem("code=proj4:+proj=utm .... +no_defs")

or from their EPSG code:

cs2 = ilwis.CoordinateSystem("code=epsg:23035")

or from a file:

cs3 = ilwis.CoordinateSystem("countries.csy")

CoordinateSystem objects themselves can do simple comparisons like:

cs1 == cs2
cs2 != cs3

CoordinateSystems can also give you the WKT or proj4 definition of the system in a string.

cs1.toWKT()
cs2.toProj4()

They also have several function for Envelopes. For example you can get and set the envelope of the coordinate system or even transform an envelope to another coordinate system.

#set the envelope of cs2 to the same as the envelope of cs1
cs2.setEnvelope(cs1.envelope())

#transform the envelope of cs3 to be in the cs2 system
env = cs3.convertEnvelope(cs2, cs3.envelope())

You can do the same transformation for a single Coordinate with coord2coord. Instead of passing a whole Envelope, just pass a single coordinate as second parameter.

CoordinateSystems are used in Coverages, Geometries and Georeferences.

FeatureCoverage

A FeatureCoverage is a collection of features. Each Feature has one or more geometries and attributes. Commonly there is one Geometry per Feature; but multiple Geometries per Feature are supported, allowing e.g. for time-series. The Geometry can be set with the geometry(...) function, and accessed by the geometry() function. For defining a Feature with multiple geometries, use the optional index parameter of those functions.

To create a Feature, use the function newFeature() on an existing FeatureCoverage object. This will create it and immediately add it to the FeatureCoverage. The parameter supplied must be a Geometry. A Geometry can be taken from another FeatureCoverage, or generated by an operation, or created by supplying a WKT string. To count the features in the FeatureCoverage, use the function featureCount().

world.featureCount()
177
mygeom = ilwis.Geometry("POINT(6.0 52.0)", world.coordinateSystem())
mygeom.toWKT()
'POINT (6.0000000000000000 52.0000000000000000)'
mygeom.coordinateSystem().name()
'countries.csy'
myfeature = world.newFeature(mygeom)
world.featureCount()
178
world.store("countries.shp", "ESRI Shapefile", "gdal")

Depending on the capabilities of the output format chosen, the last line world.store() will store the coverage as one or more datasets, as not all formats support storing mixed feature types (points, lines and polygons) in the same dataset. Attributes can be set with the setAttribute() function, and accessed with the attribute() function. The Python way of accessing attributes (myfeature["<attribute_name>"]) is also supported.

# check the value of pop_est on myfeature (the second parameter of this function is the defaultValue)
myfeature.attribute("pop_est", 0)
-1e+308
# check the value of pop_est on myfeature, python way
myfeature["pop_est"]
-1e+308
# set the value of pop_est on myfeature, python way
myfeature["pop_est"] = 1000
# confirm that it is set, python way
myfeature["pop_est"]
1000
# re-check the value, functional way
myfeature.attribute("pop_est", 0)
1000
# set the value of pop_est on myfeature, functional way
myfeature.setAttribute("pop_est", 2000)
# re-check the value, python way
myfeature["pop_est"]
2000
# inspect the current geometry
myfeature.geometry().toWKT()
'POINT (6.0000000000000000 52.0000000000000000)'
# set the geometry to a new one
myfeature.geometry(ilwis.Geometry("POINT(6 53)", world.coordinateSystem()))
# re-inspect the geometry
myfeature.geometry().toWKT()
'POINT (6.0000000000000000 53.0000000000000000)'

The attribute accessed (e.g. pop_est) must exist. These functions will not automatically add new attributes to the attribute table of the FeatureCoverage object. Use the addColumn() function on the feature coverage for this:

# attempt setting a value on a non-existing attribute
myfeature["mynewvalues"] = 999
myfeature["mynewvalues"]
...
IndexError: No attribute 'mynewvalues' found
# inspect number of columns
len(world.attributeTable().columns())
64
world.addColumn("mynewvalues", "value")
# re-inspect number of columns
len(world.attributeTable().columns())
65
# check if myfeature now has a value on the new attribute
myfeature["mynewvalues"]
-1e+308
# yes it does, -1e+308 is the default for value attributes
# set it to something else
myfeature["mynewvalues"] = 999
# read it back
myfeature["mynewvalues"]
999

FeatureIterator

Use the python way to iterate over the features. Then use the attribute get/set functions to change the feature's attributes, or use geometry() to access the geometry of the feature.

for country in world:
  country["mynewvalues"] = 1
  country.geometry().toWKT()

Geometry

This object is a placeholder for the type and the geometry data (points, vertices of lines, vertices of polygons) of a feature. The function toWKT() will reveal the string representation of the data, and fromWKT() will replace the data with the geometry represented by the supplied WKT string. Use coordinateSystem() to get the coordinate system that is used, and use setCoordinateSystem() to overwrite the coordinate system with a new one. Note that this does not alter/transform the geometry; it simply receives a new property. E.g. this can be used to overwrite a wrong coordinate system, which was e.g. the result of unsuccessfully reading the metadata of a file, or of bad or incomplete metadata. To actually transform the geometry to a different coordinate system (assuming that the current coordinate system is correct), use the function transform(). This will compute new point or vertex locations. The function envelope() returns the bounding-box of the Geometry, and the functions getArea() and getLength() return the area of polygons or the length of lines/polygons. Use distance() to compute the distance from another Geometry object. There are several functions to test the relationship between the current Geometry object and another one (e.g. within(), contains(), touches() etc.), and also several functions that will compute a new Geometry object based on only the current geometry (buffer(), convexHull()) or based on the current Geometry and another one (intersection(), Union(), difference(), symDifference()).

RasterCoverage

# open a raster coverage object:
raster = ilwis.RasterCoverage("aaraster1.mpr")
# probe its size
raster.size().xsize
1800
raster.size().ysize
1380
# get the value of a specific pixel
raster.pix2value(ilwis.Pixel(600,700))
12.201
# store it in GeoTIFF format
raster.store("aaraster", "GTiff", "gdal")

Opening netCDF / HDF5 bands in the RasterCoverage

For the "container" formats netCDF and HDF5, you must specify the layer name after the filename, separated with a slash ('/') when opening the RasterCoverage

raster = ilwis.RasterCoverage("INTER_OPER_R___TAVGD___L3__20100101T000000_20100102T000000_0003.nc/lat")
raster.size().xsize
273
raster.size().ysize
319
raster.size().zsize
1
raster.pix2value(ilwis.Pixel(0,0))
50.70261001586914

Georeference

A georeference defines the relation between rows and columns in a raster map and XY-coordinates. The location of pixels in a raster map is thus defined by a georeference. It is advised that raster maps of the same area use the same georeference. A georeference uses a coordinate system which may contain projection information. Polygon, segment and point maps merely use a coordinate system.

Every georeference has metadata like its name and size. Which can easily be accessed and altered by corresponding method calls

# check the georeference of the raster
raster.geoReference().name()
'ethtemp.grf'
# load a georeference from disk
georef = ilwis.GeoReference("Cochabamba.grf")
# inspect the size of this georeference
georef.size().xsize
637
georef.size().ysize
932

You can create a new georeference purely from scratch with a string. This string needs to contain certain key-value pairs. The keys, to which you need to provide a value are:

• type
• csy (coordinate system)
• envelope (coordinates of uper left and lower right corner)
• gridsize (x and y)
• name

# create a new georeference
# the following georeference is for a 1800x1380 grid, and the provided lat/lon bounds are assigned to the corners of this grid
georef = ilwis.GeoReference("code=georef:type=corners,csy=epsg:4326,envelope=32.991677775048 14.900003906339 47.991678557359 3.400003306568,gridsize=1800 1380,cornerofcorners=yes,name=ethnew")

After you have created a new georeference you can assign it to a raster by using setGeoReference().

# assign this new georeference to the raster
# note that the gridsize of the georeference must match the gridsize of the raster, otherwise a new empty raster is created with the new gridsize
raster.setGeoReference(georef)
raster.geoReference().name()
'ethnew'
# note that this has only altered the georeference of the in-memory RasterCoverage object, not its source file on disk.
# store the raster in GeoTIFF format, with the new projection
raster.store("aaraster", "GTiff", "gdal")

PixelIterator

The PixelIterator is a really important part of IlwisObjects. Nearly every time you modify the data of a RasterCoverage the iterator is used. You can simply create it by calling the native Python function iter() and give it an existing coverage as parameter. The iterator you get back is automatically connected to all the values fo the raster and has a size according to the size of the raster.

rc = ilwis.RasterCoverage("n000302.mpr")
#get the iterator
rcIter = iter(rc)

Basic Operations

Like every iterator the PixelIterator can be moved over the the RasterCoverage by simply adding integer values. You can also move it with Python's next() function, this will move the iterator one step forward an return the value at which it was before it moved.
Another way to get the value from the momentary position is to use the float() function. If you use the int() function you will get back the linear position of the iterator and if you use .position() you will also get the position but split up in x, y and z.
If you want to move the PixelIterator to a specific pixel you can use [] and pass it the pixel you want it to move to.

rcIter.position() #start position is pixel(0,0,0)
int(rcIter) #start linear position is 0

#move the iterator 4 forwards
rcIter += 4

#get the value of the 5th pixel (the momentary position)
float(rcIter)

#get the value and move the iterator one forward
next(rcIter)

#move iterator to pixel at (324, 124, 0)
rcIter[ilwis.Pixel(345, 124, 0)]

There is one more way to get the value of a certain pixel and that is with help of the [] operator. Just give it the linear position at which the pixel is and it will give you back its value. But you can also use this operator to set a pixel to a new value, by simply using = and then assigning the new value.

#get the value at position 475
rcIter[475]
#set the value at position 475 to 934.0
rcIter[475] = 934.0

Box

Every PixelIterator has a Box, which indicates how many pixels it can access. In general the box size equals the raster size. You can use the .box() function to get the iterators box.

bo = rcIter.box()
# get the size split up in x, y, z
bo.size()
#get the linear size
bo.size().linearSize()

begin() and end()

Apart from the iter() fucntion there are two more ways to get a PixelIterator from a RasterCoverage. If you want to get a PixelIterator that is set to the beginning of your raster you can use the function .begin(). If you want an iterator at the end of your raster, e.g. if you want to loop over a coverage using a while loop and comparing iterators, yuo can use the .end() function.

rc = ilwis.RasterCoverage("n000302.mpr")

itBegin = rc.begin()
itEnd = rc.end()

#the following loop will print all values of the iterator
while(itBegin != itEnd):
     print(next(itBegin))

You can find more on the PixelIterator in this part of the wiki.

Envelope, Coordinate, Pixel (PixelD), Size

Numpy

The first brief implementation for numpy array support uses the fromiter() factory method to create a numpy array using our PixelIterator via the standard Python iter interface. Internally ILWIS uses 'double' (8 byte floating point) representation per pixel, which maps to the np.float type.

import numpy as np
np_linear_array = np.fromiter(iter(<ilwis_raster>), np.float, <ilwis_raster>.size().linearSize())

The resulting numpy array is only a linear copy of the 2D (or even 3D) raster. The dimensions can be restored using the reshape method like this (2D):

np_linear_array.reshape( ( <ilwis_raster>.size().ysize, <ilwis_raster>.size().xsize ) )

Operations

Raster operations

The basic mathematic operations (+, -, *, /) can be executed in the Python way:

raster = ilwis.RasterCoverage("PR_2003_allyear.bsq")
raster2 = raster / 2
# this will divide all pixel values of all bands of "raster" by 2, and assign the result to "raster2"

More complicated mathematic operations (sin, cos, log ...), and all other raster operations (resample ...) are accessible through the Engine.do construct.

raster = ilwis.RasterCoverage("PR_2003_allyear.bsq")
raster2 = ilwis.Engine.do("sin", raster)
# this will compute the sinus of all pixel values of all bands of "raster", and assign the result to "raster2"

Resampling (example)

The resample operation will transform the geometry of a raster to another geometry (in other software this operation might be called "transform" or "warp")

# open a RasterCoverage
raster = ilwis.RasterCoverage("PR_2003_allyear.bsq")
# the given RasterCoverage is not georeferenced; create a suitable georeference and assign it
georef = ilwis.GeoReference("code=georef:type=corners,csy=epsg:28992,envelope=12621.630034 620583.045776 278621.630034 305583.045776,gridsize=266 315,cornerofcorners=yes,name=PR")
raster.setGeoReference(georef)
# create a new georeference; this will be the target of the `resample` operation
georef2 = ilwis.GeoReference("code=georef:type=corners,csy=epsg:28992,envelope=8548.506 621865.797 281548.506 302865.7972,gridsize=273 319,cornerofcorners=yes,name=PR_2003_warped")
raster2 = ilwis.Engine.do("resample", raster, georef2, "nearestneighbour")
# allowed values for the 3rd parameter are "nearestneighbour", "bilinear" and "bicubic"
# store raster2 to disk
raster2.store("PR_2003_warped", "GTiff", "gdal")

Cross all bands of a RasterCoverage with a PointCoverage (example)

The pointrastercrossing operation will go through all points of a PointCoverage, and for each of them, fetch the pixel value at the point's location, for all raster bands. The output is a new PointCoverage. The attribute table of the PointCoverage is the same as the original, but columns are added, one for each band of the raster.

# open a RasterCoverage
raster = ilwis.RasterCoverage("PR_2003_allyear.bsq")
# the given RasterCoverage is not georeferenced; create a suitable georeference and assign it
georef = ilwis.GeoReference("code=georef:type=corners,csy=epsg:28992,envelope=12621.630034 620583.045776 278621.630034 305583.045776,gridsize=266 315,cornerofcorners=yes,name=PR")
raster.setGeoReference(georef)
# open a PointCoverage (a pointmap); note that this must be a pointmap, thus not a polygon or segment map
pointmap = ilwis.FeatureCoverage("cities.mpp")
# perform the pointrastercrossing operation
pointmap_crossed = ilwis.Engine.do("pointrastercrossing", pointmap, raster)
# inspect the nr of columns of the original and the result
pointmap.attributeTable().columnCount()
30
pointmap_crossed.attributeTable().columnCount()
395
# the raster had 365 bands; those correspond to the added columns in pointmap_crossed.attributeTable()
raster.size().zsize
365

Empty Rasters

Empty raster are useful when generating continuous geographic data from vector data, e. g. euclidean distance around a geographical point. The code below defines a function that allow to create an empty raster with 1 band. It uses as arguments the corners of the raster, the pixel size and the spatial reference system:

# importing relevant modules
import ilwis
import math

def createEmptyRaster(ULc, LRc, pixel_size, srs):
    """
    :param ULc: a 2D tuple containing coordinates of Upper Left corner
    :param LRc: a 2D tuple containing coordinates of Lower Right corner
    :param pixel_size: size of the pixel in units of spatial reference system
    :param srs: EPSG of Spatial Reference System
    :return: an empty raster with 1 band and pixel value set to no data.
    """
    # compute raster dimensions based on resolution. Floating values are rounded UP
    width = math.ceil(math.fabs( (ULc[0] - LRc[0]) / pixel_size))
    height = math.ceil(math.fabs( (ULc[1] - LRc[1]) / pixel_size))

    # compose georeference for the raster
    georef =ilwis.GeoReference("code=georef:type=corners,csy=epsg:" + str(srs) + ",envelope=" + str(ULc[0]) + 
    " " + str(ULc[1]) + " " + str(LRc[0]) + " " + str(LRc[1]) + ", gridsize=" + str(width) + " " 
    + str(height) + ",cornerofcorners=yes,name=NoNameSet")
    # NoNameSet will be returned when calling raster.geoRegference().name()

    # Data definition for values in the band
    nr = ilwis.NumericRange(0.0, 1000000.0, 0.0)# MinValue = 0.0,  MaxValue = 1 M, precision 14 decimal places 
    #Create empty Numeric Domain
    numdom = ilwis.NumericDomain()
    #Set range of new domain
    numdom.setRange(nr)
    df = ilwis.DataDefinition(numdom, nr)

    # instance of raster
    raster = ilwis.RasterCoverage() 
    raster.setGeoReference(georef) # set size, and resolution of raster
    
    #set data definition
    raster.setDataDef(df)

    return raster

Usage: create an empty raster anchored at (231700, 485400) and (231500, 485100), resolution of 1.0 and projection RD New (EPSG: 28992). Units are defined by the spatial reference system used, in this case meters.

myEmptyRaster = createEmptyRaster(
        (231700, 485400), #ULc
        (231500, 485100), #LRc
        1.0, #resolution (pixel size)
        28992) # EPSG

Gridding (example)

Data Definition

Every coverage has a DataDefinition. With the help of a Domain it defines the type and range of data that is conatined within a raster. It is important to set it to the right range, else IlwisObjects may not be able to process the data contained in the raster and treats it as undefined values. DataDefinitions are created by calling the constructor DataDefinition and passing it a Domain and optionally a Range. You can find more on the creating of DataDefinitions in this part of the wiki.

 #Create range
 nr = ilwis.NumericRange(0.0, 10000.0)
 #Create empty Numeric Domain
 numdom = ilwis.NumericDomain()
 #Set range of new domain
 numdom.setRange(nr)

 datdef = ilwis.DataDefinition(numdom, nr)

Column Definition

Every column in a table has a own definition for its data as well. A ColumnDefinition consists of a DataDefiniton, a name and an index number. The DataDefinition defines what data can be found in this column, the name is the general name of the column and the index indicates the at which position within the table the column is positioned.

The ColumnDefinition can either be created by passing it these three objects, or instead of a DataDefinition you can as well pass a Domain

#create a new range and domain
numRan = ilwis.NumericRange(0.0, 500.0)
numDom = ilwis.NumericDomain()
numDom.setRange(numRan)
#create a new datadefinition
datDef = ilwis.DataDefinition(numDom)

#create the columndefinition
colDef = ilwis.ColumnDefinition("MARCHc", datDef, 3)
#the same effect:
colDef2 = ilwis.ColumnDefinition("MARCHc", numDom, 3)

Domains

Domains let you specify a custom range of data from existing values or newly created ones. For this to work every domain needs a specific range. The ranges differ for every kind of domain.

Every domain, except for the ColorDomain, can have a parent and accordingly a child. You can set a domain's parent with setParent(). The relation only works under the circumstances that the range of the child domain is contained within the range of the parent domain. This way you can make quick associations and link sets of data that are similar or belong thematically together.

Every domain has a method contains() which lets you check if certain data is contained within the domain's range. It returns a string depending of the containment state. These are: cNONE, if the data isn't contained, cSELF, it is contained within the domain or cPARENT, if the data is contained within the domain's parent.

Numeric Domains

For the Numeric Domains you need to specify a Numeric Range, which consists of a start value, an end value and a resolution. The resolution defines the step size in the range. If you don't define a resolution it will automatically be set to 0, which results in the highest possible resolution (as precise as a double can get).

#Create range from 5 to 60 with steps of 5. (5,10,15,20,25,...,60)
nr = ilwis.NumericRange(5, 60, 5)
#Create empty Numeric Domain
numdom = ilwis.NumericDomain()
#Set range of new domain
numdom.setRange(nr)

Item Domains

Item Domains in general are, as the name suggests Domains, whose Ranges are made up of several items. These items can have different forms. In our case these are Interval Items, Thematic Items and Named Items.
All Item Domains have the same constructor: ItemDomain() As input it takes a range and from this range it will recognize the subtype of ItemDomain it will become. You also have the possibility to give the Domain a theme, e.g. "Lake Vegetation". Which becomes importatn if you want to assign a parent, because parent and child need to have the same theme. Item Domains also give you the possibility to add and remove items. This is simply done by .removeItem() and .addItem(). For the remove function you have to provide the name of the item you want to remove. For the add function you have to provide a new valid item, depending on the Range the Domain got.

Numeric Item Domains

Numeric Item Domains are pretty similar to Numeric Domains. They also have a maximum, minimum and resolution. But additionally you can give them a name. A good use for these Domains would e.g. be to assign vegetation zones to certain heights:

# create a new NumericItemRange and fill it with items
interrange = ilwis.NumericItemRange()
interrange.add(("Tierra Caliente", 0.0, 2500.0, 1.0))
interrange.add(("Tierra Templada", 2500.0, 6000.0, 1.0))
interrange.add(("Tierra Fria", 6000.0, 12000.0, 1.0))

#assign the range to a new ItemDomain
interdom = ilwis.ItemDomain(interrange)

# set a theme
interdom.setTheme("Mountain Vegetation Zones")

# add an item
interdom.addItem(("Tierra Helada", 12000.0, 15000.0, 1.0))

# remove an item
interdom.removeItem("Tierra Fria")

# get count (will be 3 in this case) 
interdom.count()

Named Identifier Domains

Named Identifier Domains simply contain Strings as items. The according Range is called NamedItemRange.

# create and fill NamedItemRange with strings
nr = ilwis.NamedItemRange()
nr.add("Perth")
nr.add("Darwin")

# assign Range to new Domain
namedIdDom = ilwis.ItemDomain(nr)

# set a theme
namedIdDom.setTheme("Australian Cities")

# add an item
namedIdDom.addItem("Cairns")

# remove an item
namedIdDom.removeItem("Darwin")

Thematic Domains

The items of Thematic Domains consist of three string in total. These strings are its name, code and description.

#create and fill ThematicRange
tr = ilwis.ThematicRange()

tr.add("Sea", "SE", "Area covered by the sea")
tr.add("Sand", "SA", "Area covered by sand")

# assign Range to new Domain
themdom = ilwis.ItemDomain(tr)

# set a theme
themdom.setTheme("Coast")

# add an item (item to add needs to be a tuple)
themdom.addItem(("Port", "PO", "Area covered by a port"))

# remove an item
themdom.removeItem("Sand")

Time Domains

As you might guess from the name these Domains wrap a certain period of time. As its range it uses the object TimeInterval, which simply works with the date, time and datetime types that python provides in its datetime module. As with all the Domains you can use methods like contains, isValid, setParent, setStrict. But you can also get its begin and end by simply calling begin() or end(). You can set the begin and end the same way by just passing these methods a new object from the datetime module.

import datetime as dt
#create e.g. new date objects
date1 = dt.date(2014, 2, 17)
date2 = dt.date(2016, 2, 17)

#create new TimeInterval (the Range of Time Domains)
ti = ilwis.TimeInterval(date1, date2)

#assign Range to Domain
td = ilwis.TimeDomain(ti)

# get begin and end as datetime objects
tiBeg = ti.begin()
tiEnd = ti.end()

# set new begin and end (add/substract 14 days)
newBeg = date1 + dt.timedelta(14)
newEnd = date2 - dt.timedelta(14)
ti.begin(newBeg)
ti.end(newEnd)

Color Domains

At the moment the Python API for Ilwis Objects only supports Continuous Color Ranges. Which means, that you can define two colors and set one as the beginning and the other one as the ending. In the future we also plan on supporting Color Palletes. The conatins function has a slightly different name here, it is called containsColor() and takes a Color as input.

Color

To work with Color Domains, it is important that you get to know our datatype Color.

It consists of a color model and a dictionary that holds its values. The available color models that can be set are RGBA, HSLA and CYMKA, by default it will be RGBA. The elements of the dictionary can be accessed by the getItem() method which takes a string as parameter, that is used as a key for the dictionary. What the key is depends on the color model, e.g. at RGBA the keys are "red", "green", "blue" and "alpha". It is important that, when setting the values, you alwys use double numbers and no integers (1.0 is recognized, 1 is not recognized)

#create a new color with RGBA as color model and set the values
colorRGBA = ilwis.Color(ilwis.ColorModel.cmRGBA, (220.0, 20.0, 30.0, 1.0))

#create a new color with CYMKA as color model and set the values
colorCYMKA = ilwis.Color(ilwis.ColorModel.cmCYMKA, (0.6, 0.2, 0.16, 0.6, 0.5))

#create a new color with HSLA as color model
colorHSLA = ilwis.Color(ilwis.ColorModel.cmHSLA, (100.0, 0.2, 0.16, 1.0))

Continuous Color Range

As mentioned above these ranges consist of two Colors, which mark the beginning and ending of this range. They will automatically have the same model as the colors. If you want to know the ColorModel you can call defaultColorModel() on the range. If you pass a new ColorModel to the same method you can assign it to the range, but be careful with that as it might confuse the range.

color1 = ilwis.Color(ilwis.ColorModel.cmHSLA, (100.0, 0.2, 0.16, 1.0))
color2 = ilwis.Color(ilwis.ColorModel.cmHSLA, (300.0, 0.7, 0.5, 1.0))

colRan = ilwis.ContinuousColorRange(color1, color2)

Color Domain Example

color1 = ilwis.Color(ilwis.ColorModel.cmRGBA, (220.0, 20.0, 30.0, 200.0))
color2 = ilwis.Color(ilwis.ColorModel.cmRGBA, (255.0, 80.0, 60.0, 240.0))

colRan = ilwis.ContinuousColorRange(color1, color2)

colDom = ilwis.ColorDomain()
colDom.setRange(colRan)

color3 = ilwis.Color(ilwis.ColorModel.cmRGBA, (230.0, 60.0, 50.0, 240.0))

# containement check, will return cSELF, as it is in range
colDom.containsColor(color3)

Text Domains

A very simple kind of domain is a Text Domain. It simply allows strings as values.

txt = ilwis.TextDomain()

Utilities

There are several objects that occur quite often and in almost every other class, here is an overview of the most important:

Pixel

As the name suggests it is an object that represents a pixel. Pixel objects can either be created with integer numbers or doubles. In both cases you can either declare it by giving it an x, y and z value or just the x and y values.

pix = ilwis.Pixel(3,1,0)
pixXY = ilwis.Pixel(3,1)

pixD = ilwis.PixelD(0.54, 6.89, 1.2)
pixDXY = ilwis.PixelD(0.54, 6.89)

You can get, as well as set its x, y and z values with separate methods.

pix = ilwis.Pixel(0,0,0)
# set
pix.setX(34)
pix.setY(42)
pix.setZ(5)
# get
pix.x
pix.y
pix.z

You can check if the pixel has a z value (.is3D()) and compare it to another pixel.

pix1 = ilwis.Pixel(3,1)
pix2 = ilwis.Pixel(13,41,7)

pix1.is3D() #false
pix2.is3D() #true

pix1 == pix2 #false
pix1 != pix2 #true

Coordinate

As you probably can guess, they are the IwlisObjects Python API way of representing coordinates. They are very similar to Pixels, can also be set with x, y and z values or just x and y , but always consist of dooubles as values.

cXY = ilwis.Coordinate(-1871900.47, 1874237.55)
cXYZ = ilwis.Coordinate(-1871900.47, 1874237.55, 2)

The rest of the methods are exactly like the ones from the Pixel

Size

A size object is used quite often, acutally every time you want to know the size of e.g. a RasterCoverage. It can, as well, either be initialized with integer or double values. Is has an x, y and z component, of which the z component is by default set to 1, so you don't necessarily have to set it.

sz = ilwis.Size(400, 792) # will be Size(400, 792, 1)
sz = ilwis.Size(400, 792, 10)

szD = ilwis.SizeD(200.56, 233.91)
szD = ilwis.SizeD(200.56, 233.91, 2.6)

You can also get and set the separate size components or get its linear size (xyz):

sz = ilwis.Size(0, 0, 0)
#set
sz.setXsize(230)
sz.setYsize(453)
sz.setZsize(12)

#get
sz.xsize #430
sz.ysize #453
sz.zsize #12
sz.linearSize() #1250280

You can also add, substract or compare two sizes. You can also multiply it with any double value.

sz = ilwis.Size(2, 4, 5)
sz += ilwis.Size(4, 4, 0)

sz == ilwis.Size(6, 8, 5) #true
sz != ilwis.Size(6, 8, 5) #false

sz -= ilwis.Size(1, 1, 0)
sz == ilwis.Size(5, 7, 5) #true
sz *= 2 #Size(10, 14, 10)
sz *= 1.5 #Size(15, 21, 15)

You can also check if a certain Pixel is contained within the size.

sz = ilwis.Size(2, 4, 5)
ilwis.Pixel(1, 1, 1) in sz #true
ilwis.Pixel(5, 5, 5) in sz #false
ilwis.PixelD(1.2, 1.0, 4.7) in sz #true

Box and Envelope

Boxes and envelopes in IlwisObjects always define a certain area. They are in general in a rectangular shape and cover the pixels that are between there defined corners.
A box and an envelope only differ in the fact, that a box is always working with Pixels, while the envelope always uses Coordinates. The methods they use are totally identical, you only have to take care about the input, which has to be according to its type a pixel or a coordinate.

You construct a box or envelope by passing it a value for the upper left corner (the start) and for the lower right corner(the end). But you can also create either of them from a WKT string or a Size. Note that the an envelope needs a SizeD, so it needs to have doubles as values.

bo = ilwis.Box(ilwis.Pixel(3, 4, 5), ilwis.Pixel(4, 5, 6, ))
bo = ilwis.Box("POLYGON(3.6111119 4.7 5.9,4.7 5.8 6.9)")
bo = ilwis.Box(bo.size())

env = ilwis.Envelope(ilwis.Coordinate(3.6111119, 4.7, 5.9), ilwis.Coordinate(4.7, 5.8, 6.9))
env = ilwis.Envelope("POLYGON(3.6111119 4.7 5.9,4.7 5.8 6.9)")
env = ilwis.Envelope(env.size())

You can get the minimum and the maximum corner of a box/envelope and also their size.

bo.minCorner() #Pixel(3,4,5)
bo.maxCorner() #Pixel(4,5,6)
bo.size() #Size(2, 2, 2)

env.minCorner() #Coordinate(3.6111119, 4.7, 5.9)
env.maxCorner() #Coordinate(4.7, 5.8, 6.9)
env.size() #Size(2.08889, 2.1, 2)

You can check if a Pixel/Coordinate is contained within a Box/Envelope and also if another Box/Envelope is contained in it.

bo = ilwis.Box(ilwis.Pixel(0,0,0), ilwis.Pixel(10, 10, 1))
ilwis.Pixel(2, 4, 0) in bo #true
ilwis.Pixel(34, 21, 4) in bo #false
subBox = ilwis.Box(ilwis.Pixel(3,4,0), ilwis.Pixel(8, 7, 1))
subBox in bo #true

It is also possible to compare boxes with each other and envelopes with each other.

bo == subBox #false
bo != subBox #true

Footnotes