Take me to pcigeomatics.com

PCI Geomatics Help Center

How can we help you today?

Custom Pixel Arithmetic (NDVI)

PCI Geomatics -

Overview


 

This tutorial demonstrates how to perform custom raster calculations with Geomatica’s Platform Python API and NumPy. Through this tutorial, you will learn how to read in imagery using Geomatica’s Platform Python API, access the raster as a NumPy array, perform a custom calculation with multiple channels (NumPy arrays) and write the result of the calculation to an output file.

Figure 1: Example of a custom raster calculation – Each pixel from raster A is subtracted by the associated pixel in Raster B to produce raster C

Benefits


 

In some cases (or many), a specific raster routine or calculation may not exist in Geomatica or other 3rd party libraries and therefore, requiring that the calculation be custom built. This tutorial will show developers how to use Geomatica’s Python API and numpy to perform custom raster calculations that can be added to a workflow.

Prerequisites


 

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

Data Package


 

  • Click Here to download the data required for this tutorial

Tutorial


 

To demonstrate how to build custom raster calculations using Geomatica’s python API, this tutorial will read in raw landsat-8 MS imagery (as downloaded by USGS EarthExplorer) and calculate an NDVI (normalized Difference Vegetation Index) raster from the red and near-infrared (NIR) bands.

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
import os
import numpy as np
from pci.api import gobs, datasource as ds

# set input and output files for workflow
working_dir = r'D:\tutorials\raster_calculations'  # set this path to the dir where you extracted the tutorial data
raw_l8 = os.path.join(working_dir, 'LC81950282015242LGN00\LC81950282015242LGN00_MTL.txt-MS')
ndvi_out = os.path.join(working_dir, 'landsat8_ndvi.pix')

 

In the above code block, line 2 is used to import numpy, which is an open source library for storing and performing computations on large arrays. Ultimately, the custom calculations will be performed in numpy. Line 3 imports gobs and datasource classes from Geomatica’s platform api. These classes are used to read in and write out the raster pixels, geocoding information and coordinate system.

On line 6, working_dir is a variable that contains the path to the directory where your tutorial data is stored. You will need to change this to the directory on your computer where your tutorial data is saved.

2. Open the dataset and read the dataset


 

The second step is to open the dataset, access the required channels and pixels (you can subset at this step), and access/store the input coordinate system and geocoding information in variables for use later in the script.

10
11
12
13
14
15
# open the dataset
with ds.open_dataset(raw_l8) as dataset:
    reader_l8 = ds.BasicReader(dataset, [4, 5])
    in_coord_sys = reader_l8.crs
    in_geocode = reader_l8.geocoding
    in_raster = reader_l8.read_raster(0, 0, reader_l8.width, reader_l8.height)

 

On line 11 we will use a with statement to call the ds.open_dataset(raw_l8) method, passing the raw_l8 variable as it’s argument, which holds the path to the raw Landsat-8 MS image. On line 12, you create an object called reader_l8 from the ds.BasicReader() class, which is initialized by passing the dataset variable as the first argument and a list of the channels you wish to open as the second argument. In this case, we are opening channel 4 (red) and channel 5 (NIR).

On lines 13 and 14, you will call the crs and geocoding methods via the reader_l8 object. The crs and geocode methods will return the coordinate system information and geocoding information from the raw landsat-8 image to the variables in_coord_sys and in_geocode, respectively.

Finally, on line 15 you create a raster object called in_raster by calling the method read_raster(), where the arguments define the dimensions of the raster. The first two arguments define the upper-left offset and are set to zero (0), which means no offset (start reading from the upper-leftmost pixel). The last two arguments define the width and height of the raster that should be read in. In this case, we set them equal to the width and height of the input file by setting the reader_l8.width and reader_l8.height attributes. This paragraph is really just a complicated way of saying that line 15 reads in all pixels in the channels.

3. Define NumPy arrays for each input channel


 

In this step, we will define the numpy arrays for the red and NIR channels and convert them to 32-bit float. The reason we want to convert them to 32-bit float when we perform the calculation, the values in the output array will be use the same data type as the input arrays. The NDVI results will range between +1.0 and -1.0, therefore, it is important that we retain the decimal points.

17
18
19
20
21
# convert numpy array to 32bit float
red_chan = in_raster.data[:, :, 0]
red_chan_flt = red_chan.astype(np.float32)
nir_chan = in_raster.data[:, :, 1]
nir_chan_flt = nir_chan.astype(np.float32)

 

On line 18, we access the first element in the in_raster object (in_raster.data[:, :, 0]), which is the red channel and assign it to the red_chan variable. On line 19 we convert the red_chan numpy array to a 32-bit float, using astype. Line 20 and 21 do the same thing as above, but for the NIR channel.

Note: Instead of using astype, consider using vectorize to optimize performance

4. Calculating NDVI in NumPy


 

The next step is to perform an arithmetic operation with the numpy arrays that are holding the pixels for the red channel and the NIR channel, in order to calculate NDVI.

23
24
25
26
27
# Calculate NDVI and store array in ndvi_array variable (np.errstate will ignore background values set to 0)
with np.errstate(divide='ignore', invalid='ignore'):
    ndvi_array = (nir_chan_flt - red_chan_flt) / (nir_chan_flt + red_chan_flt)

ndvi_raster = gobs.array_to_raster(ndvi_array)  # wrap array in a raster object (this doesn't make a copy of the data)

 

On line 24, we use a with statement to invoke np.errstate in order to handle division errors in the case where background pixels are set to zero (0). Many orthorectified images have background pixels set as 0 and without handling this division error, your code will exit with an error. On line 25, we perform the NDVI calculation and assign the resulting array to the ndvi_array variable.

Line 27 creates a raster object for the ndvi_array called ndvi_raster, by calling the gobs.array_to_raster() method. This formats the array so that it can be written to an output file by Geomatica’s platform API.

Note: array_to_raster() does not copy the values and is optimized for performance

5. Writing the output NDVI raster object to file


 

In this final step, we will write the results of our NDVI calculation to an output file. In this case, .pix (PCIDSK), but other formats, such as GeoTIFF are also available.

 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
39
# write the output to a file
with ds.new_dataset(ndvi_out, 'PCIDSK', '') as write_dataset:
    writer = ds.BasicWriter(write_dataset)

    # Define file dimensions and write pixel values to file
    writer.create(ndvi_raster)
    writer.write_raster(ndvi_raster)

    # Add coordinate system and geocoding to output file
    writer.crs = in_coord_sys
    writer.geocoding = in_geocode

 

On line 30, we Use the with statement, we call the ds.new_dataset() function and pass the variable ndvi_out, which contains the path to the output file we want written as the first argument, ‘PCIDSK’ as the second argument, which defines the file format and an empty string ‘’ as the third argument, which could be used to define addition file options. We then use the as statement to set all this information to the write_dataset variable.

On line 31, we create an object called writer by instantiating the ds.BasicWriter() class and passing the write_dataset variable as its only argument. On line 34 we call the create() method from the writer object to create the file on disk based on the size and datatype of the raster object ndvi_raster. On line 35, we call the write_raster method to write the pixel values to the file. And finally, on lines 38 and 39, we call the writer object’s crs and geocoding methods to write the coordinate system and geocoding information from the input file (held in in_coord_sys and in_geocode variables) to the output file.

Have more questions? Submit a request

Comments

Powered by Zendesk