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.
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:
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.
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.
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.
| 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.
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.
#load required required libraries
from ctypes import *
lib1 = cdll.LoadLibrary(r'/home/jovyan/.local/lib/libGLdispatch.so.0.0.0')
lib1 = cdll.LoadLibrary(r'/home/jovyan/.local/lib/libGLX.so.0.0.0')
lib1 = cdll.LoadLibrary(r'/home/jovyan/.local/lib/libGL.so.1.7.0')
lib1 = cdll.LoadLibrary('/opt/conda/lib/libpython3.11.so.1.0')
lib1 = cdll.LoadLibrary('/opt/conda/envs/openeo/lib/libQt5Gui.so.5')
lib1 = cdll.LoadLibrary('/opt/conda/envs/openeo/lib/libQt5Sql.so.5')
lib1 = cdll.LoadLibrary('/opt/conda/envs/openeo/lib/libQt5Concurrent.so.5')
import os
import ilwis
import numpy as np
import matplotlib.pyplot as plt
gdal connector: ilwis3 connector: Organizing data:
To check the ILWIS version activated as Python site package
ilwis.version()
'1.0 build 20230718'
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.
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.
#use the Python function os.getcwd() to retrieve the current folder used and add a sub-folder
work_dir = os.getcwd() + '/Intro_ILWISPy'
#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)
/home/jovyan/mystorage/test/notebooks/ILWISPy/Intro_ILWISPy
gdal connector: ilwis3 connector: Organizing data:
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.
ilwis.catalogItems(ilwis.it.RASTER)
('dem.tif',
'ndvi.tif',
'average_monthly_temperature.mpl',
'average_monthly_temperature_april_1.mpr',
'average_monthly_temperature_august_2.mpr',
'average_monthly_temperature_december_3.mpr',
'average_monthly_temperature_february_4.mpr',
'average_monthly_temperature_january_5.mpr',
'average_monthly_temperature_july_6.mpr',
'average_monthly_temperature_june_7.mpr',
'average_monthly_temperature_march_8.mpr',
'average_monthly_temperature_may_9.mpr',
'average_monthly_temperature_november_10.mpr',
'average_monthly_temperature_october_11.mpr',
'average_monthly_temperature_september_12.mpr',
'dem.mpr',
'mydata1.mpr')
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
ilwis.catalogItems()
('ETH_M_TOWN.shp',
'dem.tif',
'dem.tif',
'dem.tif',
'ndvi.tif',
'ndvi.tif',
'ndvi.tif',
'selected.shp',
'DESCRIPN.dom',
'ETH_admin.csy',
'Landuse2.mpa',
'Landuse2.tbt',
'LanduseTheme.dom',
'MyPoints.dom',
'MyPoints.mpp',
'MyPoints.tbt',
'average_monthly_temperature.mpl',
'average_monthly_temperature.mpl',
'average_monthly_temperature_april_1.mpr',
'average_monthly_temperature_august_2.mpr',
'average_monthly_temperature_december_3.mpr',
'average_monthly_temperature_february_4.mpr',
'average_monthly_temperature_january_5.mpr',
'average_monthly_temperature_july_6.mpr',
'average_monthly_temperature_june_7.mpr',
'average_monthly_temperature_march_8.mpr',
'average_monthly_temperature_may_9.mpr',
'average_monthly_temperature_november_10.mpr',
'average_monthly_temperature_october_11.mpr',
'average_monthly_temperature_september_12.mpr',
'dem.csy',
'dem.grf',
'dem.mpr',
'ethtemp.csy',
'ethtemp.grf',
'geology.dom',
'geology.tbt',
'mydata.grf',
'mydata1.mpr',
'temperature.dom')
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.
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.
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
Loading the raster data
raster = ilwis.RasterCoverage ('ndvi.tif')
Check the raster dimensions
print(raster.size())
Size(1152, 1152, 1)
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.
print(raster.envelope())
-4607990.248000 -4599967.494000 4600009.752000 4608032.507000
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.
coordSys = raster.coordinateSystem()
coordSys.toWKT()
'PROJCS["_ANONYMOUS_1004765",GEOCS["_ANONYMOUS_1004765",DATUM["Unknown Datum,ELLIPSOID["Clarke 1866",6378206.400000000373,294.978698213898],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER[false_easting,0],PARAMETER[false_northing,0],PARAMETER[central_meridian,20],PARAMETER[latitude_of_origin,1],PARAMETER[standard_parallel_1,-19],PARAMETER[standard_parallel_2,21],UNIT[meter,1.0]]'
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.
coordSys.toProj4()
'+proj=aea +ellps=clrk66 +lat_0=1 +lat_1=-19 +lat_2=21 +lon_0=20 +x_0=0 +y_0=0'
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
datadef = raster.datadef()
print(datadef.domain())
image
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
stats = raster.statistics(ilwis.PropertySets.pHISTOGRAM)
print(stats.histogram())
[(1.0, 1455.0), (2.0, 879.0), (3.0, 163.0), (4.0, 155.0), (5.0, 179.0), (6.0, 171.0), (7.0, 197.0), (8.0, 195.0), (9.0, 260.0), (10.0, 299.0), (11.0, 374.0), (12.0, 588.0), (13.0, 901.0), (14.0, 1289.0), (15.0, 1962.0), (16.0, 2884.0), (17.0, 4029.0), (18.0, 5625.0), (19.0, 8481.0), (20.0, 12438.0), (21.0, 16615.0), (22.0, 21082.0), (23.0, 23640.0), (24.0, 21666.0), (25.0, 17214.0), (26.0, 13385.0), (27.0, 10801.0), (28.0, 9067.0), (29.0, 7446.0), (30.0, 6345.0), (31.0, 5761.0), (32.0, 5317.0), (33.0, 5139.0), (34.0, 4725.0), (35.0, 4467.0), (36.0, 4411.0), (37.0, 4420.0), (38.0, 4331.0), (39.0, 4180.0), (40.0, 3886.0), (41.0, 3881.0), (42.0, 3755.0), (43.0, 3872.0), (44.0, 3836.0), (45.0, 3740.0), (46.0, 3823.0), (47.0, 3615.0), (48.0, 3580.0), (49.0, 3503.0), (50.0, 3237.0), (51.0, 3186.0), (52.0, 3117.0), (53.0, 2977.0), (54.0, 2943.0), (55.0, 2829.0), (56.0, 2691.0), (57.0, 2513.0), (58.0, 2488.0), (59.0, 2471.0), (60.0, 2301.0), (61.0, 2256.0), (62.0, 2145.0), (63.0, 2040.0), (64.0, 2101.0), (65.0, 2031.0), (66.0, 1943.0), (67.0, 1935.0), (68.0, 1738.0), (69.0, 1805.0), (70.0, 1755.0), (71.0, 1633.0), (72.0, 1651.0), (73.0, 1590.0), (74.0, 1648.0), (75.0, 1526.0), (76.0, 1545.0), (77.0, 1506.0), (78.0, 1502.0), (79.0, 1456.0), (80.0, 1374.0), (81.0, 1331.0), (82.0, 1418.0), (83.0, 1331.0), (84.0, 1381.0), (85.0, 1287.0), (86.0, 1290.0), (87.0, 1326.0), (88.0, 1319.0), (89.0, 1260.0), (90.0, 1252.0), (91.0, 1243.0), (92.0, 1260.0), (93.0, 1288.0), (94.0, 1226.0), (95.0, 1246.0), (96.0, 1226.0), (97.0, 1267.0), (98.0, 1251.0), (99.0, 1268.0), (100.0, 1283.0), (101.0, 1287.0), (102.0, 1224.0), (103.0, 1319.0), (104.0, 1229.0), (105.0, 1413.0), (106.0, 1256.0), (107.0, 1388.0), (108.0, 1302.0), (109.0, 1322.0), (110.0, 1374.0), (111.0, 1376.0), (112.0, 1404.0), (113.0, 1415.0), (114.0, 1459.0), (115.0, 1433.0), (116.0, 1399.0), (117.0, 1464.0), (118.0, 1354.0), (119.0, 1464.0), (120.0, 1454.0), (121.0, 1435.0), (122.0, 1533.0), (123.0, 1528.0), (124.0, 1471.0), (125.0, 1468.0), (126.0, 1521.0), (127.0, 1510.0), (128.0, 1562.0), (129.0, 1565.0), (130.0, 1531.0), (131.0, 1606.0), (132.0, 1614.0), (133.0, 1627.0), (134.0, 1700.0), (135.0, 1583.0), (136.0, 1763.0), (137.0, 1737.0), (138.0, 1793.0), (139.0, 1800.0), (140.0, 1765.0), (141.0, 1859.0), (142.0, 1948.0), (143.0, 1867.0), (144.0, 1863.0), (145.0, 1926.0), (146.0, 1917.0), (147.0, 1972.0), (148.0, 1994.0), (149.0, 2005.0), (150.0, 2004.0), (151.0, 2140.0), (152.0, 2193.0), (153.0, 2100.0), (154.0, 2150.0), (155.0, 2301.0), (156.0, 2209.0), (157.0, 2193.0), (158.0, 2194.0), (159.0, 2261.0), (160.0, 2208.0), (161.0, 2261.0), (162.0, 2316.0), (163.0, 2397.0), (164.0, 2118.0), (165.0, 2194.0), (166.0, 2208.0), (167.0, 2309.0), (168.0, 2203.0), (169.0, 2174.0), (170.0, 2194.0), (171.0, 2155.0), (172.0, 2131.0), (173.0, 1997.0), (174.0, 2011.0), (175.0, 1981.0), (176.0, 1923.0), (177.0, 1833.0), (178.0, 1810.0), (179.0, 1722.0), (180.0, 1699.0), (181.0, 1539.0), (182.0, 1491.0), (183.0, 1337.0), (184.0, 1326.0), (185.0, 1257.0), (186.0, 1153.0), (187.0, 1011.0), (188.0, 973.0), (189.0, 924.0), (190.0, 842.0), (191.0, 827.0), (192.0, 710.0), (193.0, 697.0), (194.0, 577.0), (195.0, 566.0), (196.0, 522.0), (197.0, 428.0), (198.0, 421.0), (199.0, 384.0), (200.0, 312.0), (201.0, 329.0), (202.0, 253.0), (203.0, 298.0), (204.0, 214.0), (205.0, 211.0), (206.0, 200.0), (207.0, 145.0), (208.0, 154.0), (209.0, 111.0), (210.0, 101.0), (211.0, 91.0), (212.0, 73.0), (213.0, 63.0), (214.0, 50.0), (215.0, 48.0), (216.0, 38.0), (217.0, 21.0), (218.0, 27.0), (219.0, 27.0), (220.0, 25.0), (221.0, 19.0), (222.0, 12.0), (223.0, 5.0), (224.0, 4.0), (225.0, 7.0), (226.0, 1.0), (227.0, 8.0), (228.0, 0.0), (229.0, 2.0), (230.0, 3.0), (231.0, 0.0), (232.0, 0.0), (233.0, 0.0), (234.0, 2.0), (235.0, 1.0), (236.0, 0.0), (237.0, 1.0), (238.0, 0.0), (239.0, 0.0), (240.0, 0.0), (241.0, 0.0), (242.0, 0.0), (243.0, 0.0), (244.0, 0.0), (245.0, 0.0), (246.0, 0.0), (247.0, 0.0), (248.0, 0.0), (249.0, 0.0), (250.0, 0.0), (251.0, 0.0), (252.0, 0.0), (253.0, 148.0), (-1e+308, 803014.0)]
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.
print(stats[ilwis.PropertySets.pMIN]) # minimum value on the map
print(stats[ilwis.PropertySets.pMAX]) # maximum value on the map
print(stats[ilwis.PropertySets.pMEAN]) # average value, excluding undefs - nodata
print(stats[ilwis.PropertySets.pCOUNT]) # total number of pixels including undefs - nodata
print(stats[ilwis.PropertySets.pNETTOCOUNT]) # total number of pixels excluding undefs - nodata
print(stats[ilwis.PropertySets.pSUM]) # sum of all values, excluding undefs - nodata
1.0 253.0 73.04451525501345 1327104.0 524090.0 38281900.0
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
x=[a for (a,b) in stats.histogram()][:-1] #don't plot the last tuple
y=[b for (a,b) in stats.histogram()][:-1] #don't plot the last tuple
plt.plot(x,y,label='Raster values')
plt.xlabel('Data Range')
plt.ylabel('Frequency')
plt.title('Histogram')
plt.legend()
<matplotlib.legend.Legend at 0x7f84840b2510>
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!
raster.pix2value(ilwis.Pixel(500,700)) # for column,row (x,y) respectively
174.0
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'.
pixIter = iter(raster) # creates an iterator than can move over all the pixels.
pixIterAlt = raster.begin() # same iterator
#printing all the positions of pixels with the value 222
for i in range(raster.size().linearSize()):
v = float(pixIter)
if (v == 222):
print(str(pixIter.position()))
pixIter += 1
pixel(604,598,0) pixel(586,634,0) pixel(511,705,0) pixel(512,712,0) pixel(754,724,0) pixel(793,841,0) pixel(900,853,0) pixel(897,882,0) pixel(894,886,0) pixel(895,887,0) pixel(738,989,0) pixel(735,993,0)
Loading an ILWIS raster file (extension *.mpr) containing elevation information and retrieving the meta data and descriptive statistics
raster_dem = ilwis.RasterCoverage ('dem.mpr')
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
def descriptive_statistics(raster_input):
print('Raster size: ',raster_input.size())
print()
print('Map extent: ',raster_input.envelope())
print()
coordSys_input = raster_input.coordinateSystem()
print('Coordinate system: ',coordSys_input.toWKT())
print()
print('Proj4: ',coordSys_input.toProj4())
print()
datadef_input = raster_input.datadef()
print('Data type: ', datadef_input.domain())
stats_input = raster_input.statistics(ilwis.PropertySets.pHISTOGRAM)
print()
print('Minimum: ',stats_input[ilwis.PropertySets.pMIN])
print('Maximum: ',stats_input[ilwis.PropertySets.pMAX])
print('Mean: ',stats_input[ilwis.PropertySets.pMEAN])
print('No pixels - also nodata: ',stats_input[ilwis.PropertySets.pCOUNT])
print('No pixels - no nodata: ',stats_input[ilwis.PropertySets.pNETTOCOUNT])
print('Sum all: ',stats_input[ilwis.PropertySets.pSUM])
print()
x=[a for (a,b) in stats_input.histogram()][:-1]
y=[b for (a,b) in stats_input.histogram()][:-1]
plt.plot(x,y,label='Raster Map values')
plt.xlabel('Data Range')
plt.ylabel('Data Frequency')
plt.title('Raster Histogram')
plt.legend()
To execute the definition the required raster input layer has to be selected
#execute the definition using the elevation model
descriptive_statistics(raster_dem)
Raster size: Size(807, 490, 1) Map extent: 68.500000 23.850000 108.850000 48.350000 Coordinate system: PROJCS["dem.csy",GEOCS["dem.csy",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]] Proj4: +proj=longlat +a=6378137.000000000000 +b=6356752.314245179296 Data type: value Minimum: -152.0959930419922 Maximum: 6760.35205078125 Mean: 2086.7568968194337 No pixels - also nodata: 395430.0 No pixels - no nodata: 395429.0 Sum all: 825164192.9524119
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
# load once more the original DEM statistics
stats_dem = raster_dem.statistics(ilwis.PropertySets.pHISTOGRAM)
# transform the ilwis 2 dimensional raster into a 1 dimensional listing
dem_to_numpy1D = np.fromiter(iter(raster_dem), np.float64, raster_dem.size().linearSize())
# get length of the 1 dimensional listing
print('length 1D:', len(dem_to_numpy1D))
#confirm length of listing using the number of cols and rows of the raster
print('No Pixels of Raster:', (raster_dem.size().ysize * raster_dem.size().xsize))
# transform the 1 dimensional listing back into a 2 dimensional python numpy array
dem_to_numpy2D = dem_to_numpy1D.reshape((raster_dem.size().ysize, raster_dem.size().xsize))
# get the shape of the resulting 2D numpy array
print(np.shape(dem_to_numpy2D))
# extract minimum and maximum values for vizualization thresholds from the previous calculated descriptive statistics
MinValue = stats_dem[ilwis.PropertySets.pMIN]
MaxValue = stats_dem[ilwis.PropertySets.pMAX]
# create a plot of the numpy array using MatPlotLib - Imshow
fig = plt.figure(figsize =(10, 7))
plt.imshow(dem_to_numpy2D, interpolation='none', vmin=(MinValue), vmax=(MaxValue), cmap='jet')
plt.axis('on')
plt.colorbar(shrink=0.6)
plt.title('Elevation Model over the Himalayas - Central Asia')
length 1D: 395430 No Pixels of Raster: 395430 (490, 807)
Text(0.5, 1.0, 'Elevation Model over the Himalayas - Central Asia')
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'.
multibandRaster = ilwis.RasterCoverage('average_monthly_temperature.mpl')
print(multibandRaster.size())
Size(1800, 1380, 12)
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.
multibandRaster.pix2value(ilwis.Pixel(629,720,3))
19.5
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).
pixIter = iter(multibandRaster)
#printing all the positions of pixels with the value 2.40
n = multibandRaster.size().linearSize()
for i in range(n):
v = float(pixIter)
if ( v == 2.40):
print(str(pixIter.position()))
pixIter += 1
pixel(645,200,10) pixel(642,203,10) pixel(640,192,11) pixel(645,199,11)
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.
l = 720
c = 629
z = (multibandRaster.size().zsize)
print(z)
stack_value = []
for n in range(0,(z)):
point_value = (multibandRaster.pix2value(ilwis.Pixel(c,l,n)))
stack_value.append(point_value)
print('Values extracted for selected location:', stack_value)
12 Values extracted for selected location: [17.9, 18.9, 19.7, 19.5, 19.2, 17.9, 16.7, 16.7, 17.0, 17.5, 16.9, 17.2]
Create a simple plot showing the Monthly Temperature over the pixel location selected
plt.plot(stack_value, label='pixel value')
plt.ylim([14, 22])
# change X axis tick marks, as python index starts from 0
plt.xticks(np.arange(len(stack_value)), np.arange(1, len(stack_value)+1))
plt.xlabel('Month number')
plt.ylabel('Temperature (degree C)')
plt.title('Monthly temperature extracted for a given pixel')
Text(0.5, 1.0, 'Monthly temperature extracted for a given pixel')
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.
# envelope = min X, max Y, max X, min Y, here 32.20 -21.10 35.10 -18.40 (in degree.decimals for longitude and latitude)
# gridsize = no columns, no rows
# note name of georeference - here the name of output file created further below is used
grf = ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=32.20 -21.10 35.10 -18.40, gridsize=29 27, cornerofcorners=yes, name=mydata')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1000, 1))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(29,27))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)
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
array1 = np.empty(29 * 27, dtype = np.int32) #np.float64) #29 columns and 27 lines and data type as a float 64
for i in range(len(array1)):
array1[i] = i * 1
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'
rcNew.array2raster(array1)
Get dimensions of this newly created ILWIS array
print(rcNew.size())
print(rcNew.size().zsize) # number of layers
print(rcNew.pix2value(ilwis.Pixel(1,0))) # check data value for column 1,row 0 (x,y) respectively
Size(29, 27, 1) 1 1.0
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
ilwis_raster2np = np.fromiter(iter(rcNew), np.int32, rcNew.size().linearSize())
#now we overwrite the initial variable created
ilwis_raster2np = ilwis_raster2np.reshape((rcNew.size().zsize, rcNew.size().ysize, rcNew.size().xsize))
Note the order of Z,Y, X. Note that here Z is only 1 layer, same could be applicable to a multi image stack
print(ilwis_raster2np.shape)
(1, 27, 29)
To display the numpy array use can be made of the matplotlib function imshow()
plt.imshow(ilwis_raster2np[0], interpolation='none', vmin=0, vmax=850, cmap='jet') #including the layer number, note offset = 0
plt.axis('on')
plt.colorbar(shrink=0.8)
<matplotlib.colorbar.Colorbar at 0x7f8480fa41d0>
Check value of last pixel just to show that data has been properly processed, do you agree with the value retrieved?
print(rcNew.pix2value(ilwis.Pixel(28,26))) # ilwis: cols,rows, offset = 0, data type is float
print(ilwis_raster2np[0,26,28]) # numpy: z, y, x, offset = 0, data type is integer
782.0 782
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.
max_value = rcNew.size().xsize * rcNew.size().ysize - 1 #offset = 0
rcNew = ilwis.do('setvaluerange', rcNew, 0, (max_value), 1) # output as integer values and not as float
rcNew.store('mydata', 'map', 'ilwis3')
print(max_value)
782
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
rcNew.store('mydata1.mpr')
setvaluerange: in 0.000405908 seconds
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).
rcNew.store('mydata', 'GTiff', 'gdal')
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.
fc = ilwis.FeatureCoverage('ETH_M_TOWN.shp')
This loads a feature coverage containing towns (as points) in Ethiopia. So lets examine the data
fc.featureCount()
934
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.
print(fc.featureCount(ilwis.it.POLYGON))
print(fc.featureCount(ilwis.it.POINT))
print(fc.featureTypes() == ilwis.it.POINT)
print(fc.featureTypes() == ilwis.it.LINE)
0 934 True False
As with raster the geometry of the coverage in world coordinates is determined by its coordinate system
print(fc.coordinateSystem().toWKT())
print(fc.envelope())
PROJCS["ETH_M_TOWN.shp",GEOCS["ETH_M_TOWN.shp",DATUM[" D_North_American_1927",[?],ELLIPSOID["Clarke 1866",6378206.400000000373,294.978698213898],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION[""],,UNIT[degrees,1.0]] 33.066218 3.548171 46.863397 14.528715
Features can have multiple data attributes
for i in range(fc.attributeCount()):
print(str(fc[i]))
Name: ETTOWN01CO Columnindex: 0 Domain name: integer, Domain type: ValueDomain, Range: ? ? ? Name: ETTOWN01_1 Columnindex: 1 Domain name: integer, Domain type: ValueDomain, Range: ? ? ? Name: TOWN_NAME Columnindex: 2 Domain name: Text domain, Domain type: TextDomain Name: POP Columnindex: 3 Domain name: integer, Domain type: ValueDomain, Range: ? ? ? Name: WEREDA Columnindex: 4 Domain name: Text domain, Domain type: TextDomain Name: REGION Columnindex: 5 Domain name: Text domain, Domain type: TextDomain Name: ZONES Columnindex: 6 Domain name: Text domain, Domain type: TextDomain Name: POPURB Columnindex: 7 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: POPRUR Columnindex: 8 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: POPTOT Columnindex: 9 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: RAIN Columnindex: 10 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: K_VALUE Columnindex: 11 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: SOIL Columnindex: 12 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: EROSIVITY Columnindex: 13 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: ELEVATION Columnindex: 14 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: POP_TXT Columnindex: 15 Domain name: Text domain, Domain type: TextDomain Name: RAIN_TXT Columnindex: 16 Domain name: Text domain, Domain type: TextDomain Name: ELEV_TXT Columnindex: 17 Domain name: Text domain, Domain type: TextDomain Name: EASTING Columnindex: 18 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: NORTHING Columnindex: 19 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: LAT Columnindex: 20 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: LON Columnindex: 21 Domain name: value, Domain type: ValueDomain, Range: ? ? ?
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).
vector_stats = fc.statistics('POP') # POP attribute
print(vector_stats[ilwis.PropertySets.pMIN])
print(vector_stats[ilwis.PropertySets.pMAX])
print(vector_stats[ilwis.PropertySets.pMEAN])
34.0 2084588.0 7479.686295503212
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.
fiter = iter(fc) #create iterator
fiter +=4 #move it to the 4th feature
fiter.current()['POP'] # request the value of the 'POP' attribute
4388
An alternative way of accessing the data is viewing the coverage as a indexable data array and using this model to access its data.
fiter = iter(fc)
fiter[4]['POP']
4388
Because we can iterate over the features in a feature coverage the normal loops in Python will work too
for feature in fc:
if feature['POP'] == 1614:
print(str(feature))
featureid : 0 -> ETTOWN01CO : 1|ETTOWN01_1 : 45|TOWN_NAME : Adi Awalo|POP : 1614|WEREDA : Tahtay Adiabo|REGION : TIGRAY|ZONES : WESTERN|POPURB : 12761|POPRUR : 68173|POPTOT : 80934|RAIN : 521.741|K_VALUE : 32|SOIL : 1452|EROSIVITY : 128|ELEVATION : 1396.097|POP_TXT : 1614 Inh.|RAIN_TXT : 522 [mm]|ELEV_TXT : 1396 [m] asl.|EASTING : 396885.5625|NORTHING : 1606416.875|LAT : 14.528715|LON : 38.042966
Changing attribute values is also possible by using the setAttribute method of a feature.
fiter = iter(fc)
fiter += 10
fiter.current().setAttribute('POP',1000) #or fiter[10].setAttribute(...)
print(fiter[10]['POP'])
1000
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.
attTable = fc.attributeTable()
print(attTable.recordCount())
print(attTable.columnCount())
print(attTable.cell("POP",12))
934 22 24519
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.
fiter = iter(fc)
fi = fiter[10]
geom = fi.geometry()
geom.toWKT()
'POINT (38.8836776776771345 14.1930694611648178)'
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)
geom.fromWKT('POINT (39.7 15.123)')
geom.toWKT()
'POINT (39.7000000000000028 15.1229999999999993)'
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.
fcNew = ilwis.FeatureCoverage()
fcNew.setCoordinateSystem(fc.coordinateSystem())
for feature in fc:
if feature['POP'] > 100000:
fcNew.newFeatureFrom(feature,fc.coordinateSystem())
print(fcNew.featureCount()) # how many cities above 100000?
4
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.
fcNew = ilwis.FeatureCoverage()
csy = ilwis.CoordinateSystem("code=epsg:4326") # create coordinate system
fcNew.setCoordinateSystem(csy)
fcNew.setEnvelope(ilwis.Envelope('10 30 40 70')) # define the boundaries of the coverage
fcNew.addAttribute("Name", "text") # add some attributes
fcNew.addAttribute("Price","value" )
print(fcNew.attributeCount())
2
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.
feature = fcNew.newFeature('POINT (25 60)')
feature.setAttribute('Name', 'Monkey')
feature.setAttribute('Price', 1000)
feature = fcNew.newFeature('POINT (35 40)')
feature.setAttribute('Name', 'Horse')
feature.setAttribute('Price', 500)
print(fcNew.featureCount())
fiter = iter(fcNew)
print(fiter[1]['Name']) # second feature had as name 'Horse'
2 Horse
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
fcNew.store('MyPoints', 'vectormap', 'ilwis3')
fc_pol = ilwis.FeatureCoverage('Landuse2.mpa')
for i in range(fc_pol.attributeCount()):
print(str(fc_pol[i]))
Name: OBJECTID Columnindex: 0 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: AREA_ Columnindex: 1 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: PERIMETER Columnindex: 2 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: ETLUSEG_ Columnindex: 3 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: ETLUSEG_ID Columnindex: 4 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: A Columnindex: 5 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: B Columnindex: 6 Domain name: Text domain, Domain type: TextDomain Name: C Columnindex: 7 Domain name: Text domain, Domain type: TextDomain Name: D Columnindex: 8 Domain name: Text domain, Domain type: TextDomain Name: FREQUENCY Columnindex: 9 Domain name: Text domain, Domain type: TextDomain Name: DESCRIPN Columnindex: 10 Domain name: DESCRIPN.dom, Domain type: ItemDomain, Range: thematicrange: Name: LANDUSE Columnindex: 11 Domain name: Text domain, Domain type: TextDomain Name: LUT Columnindex: 12 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: USAGETYPE Columnindex: 13 Domain name: LanduseTheme.dom, Domain type: ItemDomain, Range: thematicrange: Name: Count Columnindex: 14 Domain name: count, Domain type: ValueDomain, Range: ? ? ? Name: coverage_key Columnindex: 15 Domain name: LanduseTheme.dom, Domain type: ItemDomain, Range: thematicrange:
Retrieve the classes of the polygon map
for f in fc_pol:
print(f['coverage_key'].strip('||'))
Alpine vegetation Barren lands Bushland Crop plantations Cultivation Forest Grassland Shrubland State farms Swamps Water bodies Woodland
Fetching the attribute details for the first feature from the polygon map
fi = iter(fc_pol)
f = next(fi)
print(f)
featureid : 944 -> OBJECTID : 9931|AREA_ : 0.628|PERIMETER : 16.65|ETLUSEG_ : 12736|ETLUSEG_ID : 9548|A : 639|B : 0|C : 0|D : 0|FREQUENCY : 3|DESCRIPN : Lowland bamboo bushland|LANDUSE : 73 0 0 0|LUT : 204|USAGETYPE : Alpine vegetation|Count : 17|coverage_key : Alpine vegetation||
print(f['coverage_key'])
Alpine vegetation||
Retrieving features and storing the selected features as a new shape file
fcNew = ilwis.FeatureCoverage()
fcNew.setCoordinateSystem(fc_pol.coordinateSystem())
for feature in fc_pol:
if feature['coverage_key'].strip('||') == 'Forest' or feature['coverage_key'].strip('||') == 'Swamps':
fcNew.newFeatureFrom(feature,fc_pol.coordinateSystem())
print(fcNew.featureCount())
2
fcNew.store('selected.shp', 'ESRI Shapefile', 'gdal')
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.
tbl = ilwis.Table('geology.tbt')
for i in range(tbl.columnCount()): #printing its definition of the columns
print(str(tbl[i]))
print('Number of Records: ',tbl.recordCount())
print('Number of Columns : ',tbl.columnCount())
Name: areaclass Columnindex: 0 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: areaslide Columnindex: 1 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: areanoslide Columnindex: 2 Domain name: value, Domain type: ValueDomain, Range: ? ? ? Name: densitysl Columnindex: 3 Domain name: percentage, Domain type: ValueDomain, Range: ? ? ? Name: densitynsl Columnindex: 4 Domain name: percentage, Domain type: ValueDomain, Range: ? ? ? Name: geology Columnindex: 5 Domain name: geology.dom, Domain type: ItemDomain, Range: thematicrange: Number of Records: 8 Number of Columns : 6
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.
print(tbl.cell(5,0))
print(tbl.cell("densitynsl", 3))
rec = tbl.record(4)
print(rec[0]) #areaclass
print(rec[2]) #areanoslide
Glacial deposits 97.12 10145600.0 10132800.0
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.
recs = tbl.select('densitysl > 30 and densitynsl < 40')
print(recs)
for i in range(len(recs)):
print(str(tbl.record(recs[i])))
(5,) (3069200.0, 2392000.0, 676800.0, 77.94, 22.05, 5)
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.
tbl.setCell('densitysl', 3, 100)
tbl.setCell(3, 2, 25.99) # 4th column, densitysl
print(tbl.cell('densitysl',2))
25.99
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:
ilwis.operations()
('compare',
'iscompatible',
'junction',
'numbercondition',
'range',
'createboundsonlycsy',
'createcombinationmatrix',
'createcornersgeoreference',
'createfeaturecoverage',
'createidentifierdomain',
'createintervaldomain',
'createitemrepresentation',
'createnumericdomain',
'createpalettedomain',
'createprojectedcoordinatesystem',
'createrastercoverage',
'createrastercoverage',
'createtable',
'createthematicdomain',
'createtimeddomain',
'createtimeintervaldomain',
'createvaluerepresentation',
'addcolumn',
'addrasterband',
'groupby',
'assignment',
'connect',
'convertcolumndomain',
'copycolumn',
'iffeature',
'iffraster',
'number2string',
'property',
'rastervalue',
'saveas',
'selection',
'attributeraster',
'rasterband',
'selection',
'selection',
'serviceobject',
'setattributetable',
'setvaluerange',
'setvariable',
'property',
'tablevalue',
'tablevalue',
'coord2pixel',
'coordinate',
'pixel',
'rastersize',
'mastergeoreference',
'pixel2coord',
'resample',
'spatialrelation',
'aggregation',
'binarylogicalraster',
'binarymathfeatures',
'binarymathraster',
'binarymathtable',
'joinattributes',
'columnunaryoperation',
'mapcalc',
'mapcalc',
'mapcalc',
'mapcalc',
'mapcalc',
'mapcalc',
'tabcalc',
'tabcalc',
'tabcalc',
'tabcalc',
'tabcalc',
'tabcalc',
'sin',
'cos',
'tan',
'asin',
'acos',
'atan',
'log10',
'ln',
'exp',
'abs',
'sqrt',
'ceil',
'floor',
'sgn',
'cosh',
'sinh',
'sin',
'cos',
'tan',
'asin',
'acos',
'atan',
'log10',
'ln',
'exp',
'abs',
'sqrt',
'ceil',
'floor',
'sgn',
'cosh',
'sinh',
'stringfind',
'stringsub',
'stringreplace',
'addstrings',
'testoperation',
'text2output',
'setworkingcatalog',
'buffer',
'line2polygon',
'pointrastercrossing',
'polygon2line',
'raster2point',
'transformcoordinates',
'gridding',
'MapFillSinks',
'MapFlowDirection',
'MapFlowAccumulation',
'normalizerelativedem',
'relativedem',
'script',
'aggregaterasterstatistics',
'aggregaterasterstatisticsbysubset',
'areanumbering',
'clusterraster',
'crossrasterclassification',
'cross',
'crosswithraster',
'densifyraster',
'distanceraster',
'line2raster',
'linearrasterfilter',
'lookupindex',
'mannkendallsignificancetest',
'mirrorrotateraster',
'movingaverage',
'percentilegroups',
'percentilegroupscolumn',
'point2raster',
'polygon2raster',
'rankorderrasterfilter',
'classification',
'sliceraster',
'raster2point',
'setrepresentation',
'stackminmaxpick',
'timesat',
'zonalstatisticstable',
'zonalstatisticsraster',
'aggregateraster',
'linearstretch')
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
operations = ilwis.operations()
operations = list(operations)
operations = np.unique(operations)
maxlen = max([len(operation) for operation in operations])
operations = ('\n'.join('{:3d}: {} {} {}'.format(k[0], k[1], ' ' * (2 + maxlen - len(k[1])), ilwis.operationMetaData(k[1])) for k in enumerate(operations)))
print(operations)
0: MapFillSinks MapFillSinks(inputraster,method=fill|cut) 1: MapFlowAccumulation MapFlowAccumulation(inputraster) 2: MapFlowDirection MapFlowDirection(inputraster,method=slope|height,useparalleldrainagecorrection=yes|no) 3: abs abs(raster|number); ilwis://operations/abs(table,columnname) 4: acos acos(raster|number); ilwis://operations/acos(table,columnname) 5: addcolumn addcolumn(table,columnname, domain)) 6: addrasterband addrasterband(inputraster, band) 7: addstrings addstrings(text1, text2) 8: aggregateraster aggregateraster(inputgridcoverage,aggregationmethod=!Avg|Max|Med|Min|Prd|Sum, groupsize,changegeometry) 9: aggregaterasterstatistics aggregaterasterstatistics(inputraster,statisticalmarker=mean|variance|standarddev|totalsumsquares|skew|kurtosis|max|min|maxindex|minindex|median|sum) 10: aggregaterasterstatisticsbysubset aggregaterasterstatisticsbysubset(inputraster,statisticalmarker=mean|variance|standarddev|totalsumsquares|skew|kurtosis|max|min|maxindex|minindex|median|sum, startBand, endBand) 11: aggregation aggregation(input-table, column-name|number,column-name|number) 12: areanumbering areanumbering(inputgridcoverage,connectivity=!4|8) 13: asin asin(raster|number); ilwis://operations/asin(table,columnname) 14: assignment assignment(thing) 15: atan atan(raster|number); ilwis://operations/atan(table,columnname) 16: attributeraster attributeraster(coverage,attributecolumn) 17: binarylogicalraster binarylogicalraster(gridcoverage1,gridcoverage2|number|boolean,logicaloperation=and|or|xor|less|lesseq|neq|eq|greater|greatereq) 18: binarymathfeatures binarymathfeatures(featurecoverage1,featurescoverage2,featureoperation=!add|subtract) 19: binarymathraster binarymathraster(gridcoverage|number,gridcoverage|number,rasteroperation=!add|subtract|divide|times|mod|power) 20: binarymathtable binarymathtable(input-table, !column-name|number,!column-name|number, !add|subtract|divide|times|mod) 21: buffer buffer(coverage,distance, nQuadrantSegments, endCapStyle=!round|butt|square) 22: ceil ceil(raster|number); ilwis://operations/ceil(table,columnname) 23: classification classifcation(type,multibandraster,sampleraster,widen-factor|threshold-distance) 24: clusterraster clusterraster(inputraster,numberOfCluster[,createAttributeTable]] 25: columnunaryoperation columnunaryoperation(input-table, column-name|number, Avg|Max|Med|Min|Pred|Std|Sum[,output-column-name]) 26: compare compare(firstvalue,operator=!equals|notequals|greatherthan|lessthan|greaterorequal|lessorequal,secondvalue) 27: connect connect(url) 28: convertcolumndomain convertcolumndomain(inputtable,columnname,targetdomaintype=!identifier|thematic|time|float|integer|color, domain-name) 29: coord2pixel coord2pixel(datatype=rastercoverage|georef,Coordinate) 30: coordinate coordinate(x,y[,z]) 31: copycolumn copycolumn(sourcetable,sourcecolumn, targettable, targetcolumn, sourceKey, targetKey)) 32: cos cos(raster|number); ilwis://operations/cos(table,columnname) 33: cosh cosh(raster|number); ilwis://operations/cosh(table,columnname) 34: createboundsonlycsy createcoundsonlycsy(minx, miny, maxx, maxy) 35: createcombinationmatrix createcombinationmatrix(domainxaxis,domainyaxis,domaincombos, xaxisvalues,yaxisvalues,combovalues) 36: createcornersgeoreference createcornersgeoreference(minx, miny, maxx, maxy, pixelsize,coordinatesystem,centerofpixels,description) 37: createfeaturecoverage createfeaturecoverage(coordinatesystem, coveragelist,stackdomain,stackdefinition,mergelayers 38: createidentifierdomain createidentifierdomain(itemdefintion, strict,description[,parentdomain]) 39: createintervaldomain createintervaldomain(intervaldefinition, resolution, strict,description[,parentdomain]) 40: createitemrepresentation createitemrepresentation(basedomain, items) 41: createnumericdomain createnumericdomain(min,max,resolution, strict,description[,parentdomain])) 42: createpalettedomain createpalettedomain(itemdefintion, strict,description[,parentdomain]) 43: createprojectedcoordinatesystem createprojectedcoordinatesystem(projectionname,projection parameters,ellipsoid[,datumshifts][,envelope]) 44: createrastercoverage createrastercoverage(georeference, domain,bands,stackdomain=positiveInteger, stack-defintion=1..,auto-resample); createrastercoverage(georeference, empty=!true|false) 45: createtable createtable(columns)) 46: createthematicdomain createthematicdomain(itemdefintion, strict,description[,parentdomain]) 47: createtimeddomain createtimedomain(min,max,resolution, description) 48: createtimeintervaldomain createtimeintervaldomain(intervaldefinition, resolution,description) 49: createvaluerepresentation createvaluerepresentation(basedomain, items, relative, steps) 50: cross cross(raster1, raster2, undefhandling=!ignoreundef|ignoreundef1 | ignoreundef2 | dontcare) 51: crossrasterclassification crossrasterclassification(inputraster1,inputraster2, combinationmatrix) 52: crosswithraster crosswithraster(raster1, raster2, undefhandling=!ignoreundef|ignoreundef1 | ignoreundef2 | dontcare) 53: densifyraster densifyraster(raster, factor, interpolation=!nearestneighbour | bilinear | bicubic) 54: distanceraster distanceraster(raster, weightraster) 55: exp exp(raster|number); ilwis://operations/exp(table,columnname) 56: floor floor(raster|number); ilwis://operations/floor(table,columnname) 57: gridding gridding(coordinatesystem,top-coordinate,x-cell-size, y-cell-size, horizontal-cells, vertical-cells) 58: groupby groupby(inputtable,aggregationColumn, aggregationMethod=!sum|average|maximum|minimum|predominant|last|first) 59: iffeature iffeature(featurecoverage,outputchoicetrue, outputchoicefalse) 60: iffraster iffraster(rastercoverage,outputchoicetrue,outputchoicefalse) 61: iscompatible iscompatible(firstvalue,secondvalue,georeference|coordinatesystem|domain|raster|featurecoverage|table) 62: joinattributes joinattributes(base-table,column-name,ret-column, input-table, column-name, ret-columns) 63: junction junction(FirstValue, SecondValue, Type=!RasterCoverage|FeatureCoverage|Number|String|Table|Boolean) 64: line2polygon line2polygon(inputfeatures) 65: line2raster line2raster(inputlinecoverage,targetgeoref ) 66: linearrasterfilter linearrasterfilter(raster1, linearfiltername) 67: linearstretch linearstretch(raster,percentage) | linearstretch(raster,percentage,minout,maxout) | linearstretch(raster,minin,maxin) | linearstretch(raster,minin,maxin,minout,maxout) 68: ln ln(raster|number); ilwis://operations/ln(table,columnname) 69: log10 log10(raster|number); ilwis://operations/log10(table,columnname) 70: lookupindex lookupindex(raster, testraster, operator=!smaller|smallerequal|equal|notequal|greater|greaterequal) 71: mannkendallsignificancetest mannKendallSignificanceTest(raster,domain) 72: mapcalc mapcalc(expression,gridcoverage|number); mapcalc(expression,gridcoverage|number,gridcoverage|number); mapcalc(expression,gridcoverage|number,gridcoverage|number,gridcoverage|number); mapcalc(expression,gridcoverage|number,gridcoverage|number,gridcoverage|number,gridcoverage|number); mapcalc(expression,gridcoverage|number,gridcoverage|number,gridcoverage|number,gridcoverage|number,gridcoverage|number); mapcalc(expression,gridcoverage|number,gridcoverage|number,gridcoverage|number,gridcoverage|number,gridcoverage|number,gridcoverage|number) 73: mastergeoreference mastergeoref(targetgeoref) 74: mirrorrotateraster mirrorrotateraster(inputraster,mirrortype=mirrhor | mirrvert | mirrdiag | transpose | rotate90 | rotate180 | rotate270) 75: movingaverage movingaverage(inputpointmap,attribute,invDist | linear,exp,limDist,georef| xsize[,ysize]) 76: normalizerelativedem normalizerelativedem(inputrelativedem,inputcatchment) 77: number2string number2string(number) 78: numbercondition numbercondition(FirstValue, Condition=!greater than|smaller than|greater or equals|smaller or equals|equals|not equals, SecondValue) 79: percentilegroups percentilegroups(inputraster, distributionGroups, sizeControlGroup) 80: percentilegroupscolumn percentilegroupscolumn(inputraster, distributionGroups, sizeControlGroup,yes|!no) 81: pixel pixel(x,y[,z]) 82: pixel2coord pixel2coord(datatype=rastercoverage|georef,Pixel) 83: point2raster point2raster(inputpointmap,targetgeoref | xsize[,ysize]) 84: pointrastercrossing pointrastercrossing(inputpointmap, inputgridcoverage) 85: polygon2line polygon2line(inputfeatures, usesingleid=true) 86: polygon2raster polygon2raster(input-polygonmap,targetgeoref) 87: property property(ilwisobject,property=!name|xsize|ysize|zsize|rows|columns|type|valuetype|featurecount|pointcount|linecount|polygoncount|pixelcount|domain|projection|ellipsoid|url|rawurl|min|max); property(systemproperty=workingcatalog) 88: range range(number|string|band|polygon|line|point|feature|record,container) 89: rankorderrasterfilter rankorderrasterfilter(raster1, filtername [,columns,rows,index] 90: raster2point raster2point(inputraster); ? 91: rasterband rasterband(raster,stackindex) 92: rastersize gridsize(rastercoverage,dimension=xsize|ysize|zsize) 93: rastervalue rasvalue(inputgridcoverage,x,y,[,z]) 94: relativedem relativedem(inputdemcoverage,inputflowdircoverage,inputisdrain) 95: resample resample(inputgridcoverage,targetgeoref,interpolation=nearestneighbour|bilinear|!bicubic) 96: saveas saveas(inputobject,url,outputformat, outputprovider) 97: script script file|script scriptline(,scriptline)* 98: selection selection(coverage,selection-definition); selection(featurecoverage,selection-definition); selection(table,selection-definition) 99: serviceobject serviceobject(ilwisobject,aspect=!georeference|domain|coordinatesystem|projection|ellipsoid) 100: setattributetable setattributetable(inputtable, raster, primarykey) 101: setrepresentation stackminmaxpick(raster1, attribute=Pixel value, representation 102: setvaluerange setvaluerange(objectname,min,max,resolution) 103: setvariable setvariable(value) 104: setworkingcatalog setworkingcatalog(url) 105: sgn sgn(raster|number); ilwis://operations/sgn(table,columnname) 106: sin sin(raster|number); ilwis://operations/sin(table,columnname) 107: sinh sinh(raster|number); ilwis://operations/sinh(table,columnname) 108: sliceraster sliceraster(inputgridcoverage, IntervalDomain) 109: spatialrelation contains(featurecoverage,geometries,contains|covers|coveredBy|touches|intersects|disjoint|within|equals|crosses|overlaps) 110: sqrt sqrt(raster|number); ilwis://operations/sqrt(table,columnname) 111: stackminmaxpick stackminmaxpick(raster1, minimum|maximum 112: stringfind stringfind(source,searchtext,[,begin]) 113: stringreplace StringReplace(source,searchtext,replacetext) 114: stringsub stringsub(source,begin,[,end]) 115: tabcalc tabcalc(expression,inputtable,outputcolumn,yes|!no, inputcolumn); tabcalc(expression,inputtable,outputcolumn,yes|!no, inputcolumn, inputcolumn); tabcalc(expression,inputtable,outputcolumn,yes|!no, inputcolumn, inputcolumn, inputcolumn); tabcalc(expression,inputtable,outputcolumn,yes|!no, inputcolumn, inputcolumn, inputcolumn, inputcolumn); tabcalc(expression,inputtable,outputcolumn,yes|!no, inputcolumn, inputcolumn, inputcolumn, inputcolumn, inputcolumn); tabcalc(expression,inputtable,outputcolumn,yes|!no, inputcolumn, inputcolumn, inputcolumn, inputcolumn, inputcolumn, inputcolumn) 116: tablevalue tablevalue(inputtable,columnname,recordnr); tablevalue(inputtable,primarycolumnname,primarykeyvalue, valuecolumn) 117: tan tan(raster|number); ilwis://operations/tan(table,columnname) 118: testoperation testoperation(raster1, raster2) 119: text2output text2output(text[,text]) 120: timesat timesat(inputgridcoverage,iterationcount,upperenvelop,fitlastiteration,extendwindow) 121: transformcoordinates transformcoordinates(inputfeaturemap, csydefintion) 122: zonalstatisticsraster zonalstatisticsraster(raster1, raster2, table, column, method=!average|sum|max|min) 123: zonalstatisticstable zonalstatisticstable(raster1, raster2, table, column, method=!average|sum|max|min)
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).
print(ilwis.operationMetaData('resample'))
print()
print(ilwis.operationMetaData('resample', 'description'))
print(ilwis.operationMetaData('resample', 'inparameters')) # how many parameters needed
print(ilwis.operationMetaData('resample', 'outparameters')) # how many output values
print(ilwis.operationMetaData('resample', 'keyword'))
resample(inputgridcoverage,targetgeoref,interpolation=nearestneighbour|bilinear|!bicubic) translates a rastercoverage from one geometry (coordinatesystem+georeference) to another 3 1 raster,geometry,transformation,georeference,projection,coordinatesystem,reprojection
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.
print(ilwis.operationMetaData('resample','input',1,'name')) # first input parameter
print(ilwis.operationMetaData('resample','input',1,'description'))
print(ilwis.operationMetaData('resample','input',1,'type'))
print(" ")
print(ilwis.operationMetaData('resample','input',2,'name')) # second input parameter
print(ilwis.operationMetaData('resample','input',2,'description'))
print(ilwis.operationMetaData('resample','input',2,'type'))
print(" ")
print(ilwis.operationMetaData('resample','input',3,'name')) # third input parameter
print(ilwis.operationMetaData('resample','input',3,'description'))
print(ilwis.operationMetaData('resample','input',3,'type'))
print(" ")
print(ilwis.operationMetaData('resample','output',1,'name')) # output value
print(ilwis.operationMetaData('resample','output',1,'description'))
print(ilwis.operationMetaData('resample','output',1,'type'))
input rastercoverage input rastercoverage with domain any domain rastercoverage target georeference the georeference to which the input coverage will be morphed georeference Resampling method The method used to aggregate pixels from the input map in the geometry of the output map string output rastercoverage output rastercoverage with the domain of the input map rastercoverage
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'
raster = ilwis.RasterCoverage('ndvi.tif')
csyLatLon = ilwis.CoordinateSystem('code=epsg:4326')
latLonEnv = csyLatLon.convertEnvelope(raster.coordinateSystem(), raster.envelope())
grfTarget = ilwis.GeoReference(csyLatLon, latLonEnv,raster.size())
Now we can do the actual operation and fill in the parameters according to the syntax
rcNew = ilwis.do('resample', raster, grfTarget, 'bilinear')
print('\n')
print(rcNew.size())
print(rcNew.coordinateSystem())
print(rcNew.coord2value(ilwis.Coordinate(22.78, 5.034))) #querying a location in lon/lat instead of Albers.
resample: in 1.29856 seconds Size(1152, 1152, 1) LatLon WGS 84 171.50697222324135
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.
rcNew = ilwis.do('mapcalc', '@1 + 100', raster)
mapcalc: in 0.157027 seconds
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.
rc2 = ilwis.do('mapcalc', '(log10(@1) * 25 + @2)/ @2', raster, rcNew)
mapcalc: in 0.239583 seconds
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
rc3 = ilwis.do('mapcalc', 'iff(@1>40, ?, 1)', raster)
mapcalc: in 0.153219 seconds
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
multibandRaster = ilwis.RasterCoverage('average_monthly_temperature.mpl')
rc4 = ilwis.do('mapcalc','@1 * 10', multibandRaster)
mapcalc: in 4.0874 seconds
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
rcNew = ilwis.do('select',raster, 'boundingbox(300 250, 650 500)')
selection: in 0.00416295 seconds
This cuts(selects) a section of the raster and puts it in a new raster
rcNew = ilwis.do('select',raster, 'pixelvalue > 40 with: boundingbox(300 250, 650 500)')
selection: in 0.0116116 seconds
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.
rcNew = ilwis.do('select',raster, 'pixelvalue > 40 and pixelvalue < 150 with: boundingbox(300 250, 650 500)')
selection: in 0.0295599 seconds
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.
The geometries support a number of the operations as defined by OGC. e.g. The distance between two features
fc = ilwis.FeatureCoverage('ETH_M_TOWN.shp')
fiter = iter(fc)
f1 = fiter[10] # tenth feature
f2 = fiter[100] # 100th feature
g1 = f1.geometry()
g2 = f2.geometry()
print(g1.distance(g2))
#or
g3 = g1.buffer(0.5)
print(g3.getArea())
3.1132397050449203 0.7803612880645107
The operation supported are distance, getArea, getLength, buffer, convexHull, intersection, difference, symetricDifference and join. E.g.
csy = ilwis.CoordinateSystem("code=epsg:4326")
geom1 = ilwis.Geometry('Polygon((35.9 6.5, 35.9 8.6, 40.5 8.6, 40.5 6.5, 35.9 6.5))', csy)
geom2 = ilwis.Geometry('Polygon((34.7 6.0, 34.7 8.0, 40.0 8.0, 40.0 6.0, 34.7 6.0))', csy)
geom3 = geom1.intersection(geom2)
print(geom3.toWKT())
print(geom3.getLength())
POLYGON ((35.8999999999999986 6.5000000000000000, 35.8999999999999986 8.0000000000000000, 40.0000000000000000 8.0000000000000000, 40.0000000000000000 6.5000000000000000, 35.8999999999999986 6.5000000000000000)) 11.200000000000003
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.
print(geom1.overlaps(geom2))
print(geom1.overlaps(g1))
True False