PCI Geomatics -

# Overview

This tutorial will guide you through the process of creating a custom pixel-wise processing routine with Geomatica and Python. In this context, pixel-wise processing is considered to be running a routine on an pixel location and returning a result, then move on to the next pixel location and repeat the same routine, until all pixels are processed. For example, in the graphic below, the developer created a function that checks each pixel position in two images for the color of the pixel in each image. If they match, the output pixel is colored purple. Otherwise the pixel is colored yellow.

# Benefits

Many algorithms or workflows need to evaluate the value of a pixel at a specific location against itself or against pixel values from other images or channels at the same location. This tutorial will teach you how to create completely custom routines for processing every pixel in an image one at a time.

# 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
• How to create if-statements
• Some experience programming with Geomatica (recommended)
• Some experience with numpy (recommended)
• Understanding how to use a decorator (recommended)

# Tutorial

To demonstrate how to build custom pixel-wise processes, we will expand on the Raster Arithmetic tutorial. Specifically, we will calculate an NDVI (normalized Difference Vegetation Image) from a Landsat-8 image and then create a custom pixel-wise process that evaluates the NDVI value at each pixel. The pixel-wise process will evaluate whether a pixel contains high vegetation cover, moderate vegetation cover or no vegetation.

## 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'C:\pixel_wise_proc' # set this path to the dir where you extracted the tutorial data raw_l8 = os.path.join(working_dir, 'LC81950282015242LGN00\LC81950282015242LGN00_MTL.txt-MS') veg_out = os.path.join(working_dir, 'landsat8_veg_map.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``` ```# Calculate NDVI and store array in ndvi_array variable (np.errstate will ignore background values set to 0) print '\nCalculating NDVI Image...' with np.errstate(divide='ignore', invalid='ignore'): ndvi_array = (nir_chan_flt - red_chan_flt) / (nir_chan_flt + red_chan_flt) ```

On line 25, 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 26, we perform the NDVI calculation and assign the resulting array to the ndvi_array variable.

## 5. Create function with pixel-wise process to identify vegetation cover

In this step we will create a decorator for numpy’s np.vectorize function and then create the custom function that will define how we should process each pixel.

 ``` 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43``` ```@np.vectorize def find_veg(ndvi): if 0.3 < ndvi < 0.4: return 150 # Moderate vegetation cover (Yellow) elif ndvi >= 0.4: return 100 # high vegetation cover (Green) else: return 255 # no vegetation cover (White) print '\nComputing Vegetation Results...' with np.errstate(invalid='ignore'): veg_id_array = find_veg(ndvi_array) veg_id_int = veg_id_array.astype(np.uint8) veg_id_raster = gobs.array_to_raster(veg_id_int) ```

On line 29 we create a decorator with the @ character for the np.vectorize function, which will take the custom function that directly follows (find_veg) and run that function on each value in a numpy array (each pixel).

On line 30, we create the custom function find_veg(ndvi):, which takes in the ndvi _array as its only argument (array to process).

Lines 31 to 36 contain the custom rules for determining if a pixel is high vegetation cover, moderate vegetation cover or no vegetation cover.

Line 42 converts the output array, veg_id_array, to an 8bit unsigned integer. Line 43 creates a raster object for the veg_id_int array called veg_id_raster, by calling the gobs.array_to_raster() method, which 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

## 6. Writing the output vegetation map to file

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

 ``` 45 46 47 48 49 50 51 52 53 54 55``` ```# write the output to a file with ds.new_dataset(veg_out, 'PCIDSK', '') as write_dataset: writer = ds.BasicWriter(write_dataset) # Define file dimensions and write pixel values to file writer.create(veg_id_raster) writer.write_raster(veg_id_raster) # Add coordinate system and geocoding to output file writer.crs = in_coord_sys writer.geocoding = in_geocode ```

On line 46, we Use the with statement to call the ds.new_dataset() function and pass the variable veg_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 47, we create an object called writer by instantiating the ds.BasicWriter() class and passing the write_dataset variable as its only argument. On line 50 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 51, we call the write_raster method to write the pixel values to the file. And finally, on lines 54 and 55, 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.