Take me to pcigeomatics.com

PCI Geomatics Help Center

How can we help you today?

EASI modelling in Python

PCI Geomatics -

Overview


 

This tutorial will demonstrate how users can create simple models with Geomatica’s EASI syntax (legacy programming language) and run them directly from their python scripts using Geomatica’s MODEL algorithm.

Benefits


 

EASI (Engineering and Applied Science Interface) is PCI Geomatics’ legacy scripting and programming language for accessing and automating Geomatica workflows. While there are no plans to continue development with EASI, many current customers have extensive models and workflows that they may wish to access through Python, as they begin to transition completely over to Python. Furthermore, there are some concepts that may be easier to build in EASI for those who are familiar with it.

Prerequisites


 

  • Geomatica 2015 or later installed
  • Python 101 (basic understanding of python coding concepts)
    • How to access modules, classes, attributes and methods
    • How to call a function
  • Some experience programming with Geomatica (recommended)
    • In particular, EASI programming

Data Package


 

  • Click Here to download the data required for this tutorial

Tutorial


 

To demonstrate this concept, we will build an EASI model that identifies ice/snow and separately, water bodies from Landsat-8 multispectral (MS) and thermal imagery, using the SWIR 2.2µm and thermal bands. We will then extract snow/ice objects as one polygon layer and water bodies as a second polygon layer.

1. Import necessary modules & setup inputs/outputs


 

The first step is to import the Geomatica modules we will be calling in this tutorial and setting up the variables that will point to our input file(s) and output file(s).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
import os
from pci.pcimod import pcimod
from pci.model import model
from pci.dn2toa import dn2toa
from pci.atcor_t import atcor_t
from pci.iii import iii
from pci.expolras import expolras

# set input and output files for workflow
working_dir = r'C:\easi_modelling'  # set this path to the dir where you extracted the tutorial data
raw_thermal_file = os.path.join(working_dir, 'LC81950282015242LGN00', 'LC81950282015242LGN00_MTL.txt-THERMAL')
raw_ms_file = os.path.join(working_dir, 'LC81950282015242LGN00', 'LC81950282015242LGN00_MTL.txt-MS')
toa_file = os.path.join(working_dir, 'l8_toa.pix')
temper_file = os.path.join(working_dir, 'l8_temperature.pix')
ice_snow_poly = os.path.join(working_dir, 'ice_snow_poly.pix')
water_bodies = os.path.join(working_dir, 'water_poly.pix')

 

In the above code block, line 1 will import the standard python module os. Lines 2 to 7 import the different Geomatica modules that we will be working with:

  • dn2toa – Geomatica’s algorithm for converting satellite imagery from at-sensor radiance to top of atmospheric reflectance (TOA)
  • actor_t – Geomatica’s algorithm for calculating temperature values from thermal infrared satellite imagery
  • pcimod – Geomatica’s algorithm for modifying a .pix file (adding or removing image channels)
  • iii – Geomatica’s algorithm for copying raster data (imagery) from one file to another
  • model – Geomatica’s algorithm for executing EASI scripts and models in python
  • expolras – Geomatica’s algorithm for extracting polygons from raster data

On line 10 we create a variable called working_dir. Make sure to set it equal to the path where you extracted the tutorial data (root folder). On lines 11 to 16, we create a number of variables that hold paths to various input files and output files, which are required by this workflow.

The final results from this tutorial will be written to the files in the ice_snow_poly and water_bodies variables on lines 15 and 16, respectively.

2. Compute TOA (for SWIR 2.2µm) and temperature (for TIR) images


 

Raw satellite imagery is often distributed as scaled at-sensor radiance values, whose pixel values can range significantly from one image to the next. If we create a model that evaluates a feature based on these pixel values, it would likely only work for this specific image. By converting the SWIR image to TOA, we can normalize the values to a range between 0 and 100%. The same concept is true by computing the temperature values from the TIR channel. This will help make our model more robust and reusable on other images.

18
19
20
# Calculate TOA and temperature for Landsat-8 MS and THERMAL, respectively
dn2toa(fili=raw_ms_file, filo=toa_file)
atcor_t(fili=raw_thermal_file, filo=temper_file)

 

On line 19 we call the Geomatica algorithm dn2toa(), which calculates the TOA image. The following keyword arguments are passed to the algorithm:

  • fili=raw_ms_file – sets the input file that the TOA computation should be run on
  • filo=toa_file – sets the output file that the TOA results will be written to

On line 20 we call the Geomatica algorithm actor_t(), which calculates temperature values from the TIR image. The following keyword arguments are passed to the algorithm:

  • fili=raw_thermal_file – sets the input TIR file that the temperature computations should be run on
  • filo=temper_file – sets the output file that the temperature values will be written to

 

3. Transfer TOA SWIR channel to Temperature file


 

The model() function, which is used to execute EASI code, requires that all channels used in the model are in the same .pix file. Therefore, we are required to copy the SWIR 2.2µm TOA image channel to the temperature file.

22
23
24
# Transfer temperature channel to toa file
pcimod(file=temper_file, pciop='add', pcival=[1, 0, 0, 1, 0, 0])
iii(fili=toa_file, filo=temper_file, dbic=[7], dboc=[3])

 

On line 23 we call the Geomatica algorithm pcimod(), which adds empty channels to temper_file, which is required before we can copy the SWIR 2.2µm TOA channel to temper_file. The following keyword arguments are passed to the algorithm:

  • file=temper_file – sets the file that the channels should be added to or removed from
  • pciop=’add’ – tells the algorithm that we wish to add channels
  • pcival=[1, 0, 0, 1, 0, 0] – creates an 8bit and 32bit channel

On line 24 we call the Geomatica algorithm iii() which copies an image channel from one file to another. The following keyword arguments are passed to the algorithm:

  • fili=toa_file – sets the file that we want to copy the channel from
  • filo=temper_file – sets the destination file that will receive the copied channel
  • dbic=[7] – sets the input channel number to copy (SWIR 2.2)
  • dboc=[3] – sets the output channel that will receive the channel (the 32bit channel)

 

4. Create modelling conditions in EASI syntax and run model algorithm


 

The next step is to create the conditionals for identifying ice/snow and separately, water with EASI syntax, assign it to a string variable and then run it with the Geomatica algorithm model().

 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
# Create EASI statement to identify snow and ice and run through MODEL algorithm
# Ice/snow = 5 and water/ice water = 10
my_model = '''if(%1 <> 0) and (%1 <= 5) and (%3 < 3)then
    %2=5;
elseif(%1 <> 0) and (%1 >= 20) and (%3 < 0.75) then
    %2=10;
elseif(%1 <> 0) and (%1 >= 5) and (%1 <= 20) and (%3 < 0.75) then
    %2=10;
endif;'''

# Run the MODEL algorithm
model(file=temper_file, source=my_model, undefval=[])

 

On line 28 we create a variable called my_model and set it equal to a docstring with three consecutive apostrophes, which allows us to create a string that will span multiple lines. We then put the EASI code inside the docstring (lines 28 to 34).

EASI Programming Basics: The percentage sign directly followed by a value refers to a channel (i.e. %1 is channel 1). Two percentage signs directly followed by a value refers to a segment number (i.e. %%2 is segment 2). Visit the EASI help page for more information about programming with EASI.

On line 37, we call the model() algorithm, which will be used to execute the EASI code held in my_model variable. The following keyword arguments are passed to the algorithm:

  • file=temper_file – tells model() which file to run the EASI script on
  • source=my_model – point to the file or string with the EASI code
  • undefval=[] – sets value for undefined operation (in this case no_data_value)

 

5. Extract Polygons for ice/snow and water objects, separately


 

Here we will run the Geomatica algorithm expolras() twice, the first time to extract ice/snow polygons and the second time to extract water bodies.

39
40
41
42
43
# Extract snow and ice as polygon layer
expolras(fili=temper_file, filo=ice_snow_poly, dbic=[2], tval=[5, 5], areaval=[3])

# Extract water bodies as polygon layer (including ice covered water bodies)
expolras(fili=temper_file, filo=water_bodies, dbic=[2], tval=[10, 10], areaval=[1500])

 

On line 40, we call the expolras() algorithm to extract polygons for ice and snow objects (glaciers). The following keyword arguments are passed to the algorithm:

  • fili=temper_file – Input file for extracting polygons
  • filo=ice_snow_poly – output file with the polygon layer
  • dbic=[2] – input channel for extracting polygons
  • tval=[5, 5] – value we wish to extract
  • areaval=[3] – minimum area set to 3 pixels

On line 43, we call the expolras() algorithm again and run it on the same dataset, but this time set tval=10, which is the value we used to identify water bodies. We also set areaval=1500

Have more questions? Submit a request

Comments

Powered by Zendesk