Take me to pcigeomatics.com

PCI Geomatics Help Center

How can we help you today?

Change Detection and Area Reporting

PCI Geomatics -

Overview

This tutorial outlines how to detect change between two images and then create an area report stating the percentage and total area of change.

Benefits

Change detection is useful in numerous circumstances in which you may want to analyze change over time, such as:

  • Deforestation
  • Storm damage
  • Forest-fire damage
  • Flooding
  • Urban sprawl

Prerequisites

  • Geomatica 2015 or later installed
  • Python 101 (basic understanding of python coding concepts)
  • Some experience programming with Geomatica (recommended)

Data Package

  • Click Here to download the data required for this tutorial

Tutorial

The data that will be used are two Landsat-5 TM images (1990 and 1999) of the Prince George area in Northern British Columbia, Canada. This area is well known for the logging activities that occur here. We are specifically interested in looking at areas of deforestation between these two times. The following workflow will be followed to create the change detection image and area report.

 

1990 Landsat 5 Image

1999 Landsat 5 Image

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 and output directories and files.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
from pci.chdetop import chdetop
from pci.expolras import expolras
from pci.areareport import areareport
from pci.map import map as bit_encode
from pci.pcimod import pcimod
from pci.poly2bit import poly2bit
from pci.nspio import Report, enableDefaultReport
import os

working_dir= r'I:\Data\Tutorial\Python\ChangeDetection'
input_dir = os.path.join(working_dir, "Input")
change_image= os.path.join(working_dir, "change_image.pix")

 

Lines 1 to 8 import the six required algorithms, as outlined above. You will notice that on line 4 the map algorithm is imported as bit_encode. This was done because map() is a built-in python function. It is best practice to import this Geomatica algorithm using a different name as to not override the python function. On line 10, working_dir is a variable that contains the path to the directory where your change image will reside. You will also need to set the input_dir, where the two Landsat images are stored. Finally, the change_image is set on line 12. This will be the pix file where you will save the outputs from each of the algorithms.

2. Create change detection raster and extract change polygons

The first processing step will be to create a raster which shows the areas of change and then extract those change areas as polygons. The first algorithm you will run is chdetop, which is used to run change detection on optical imagery.

 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
28
# Create change detection raster
chdetop(fili=os.path.join(input_dir,"4922_etm_12sep1999_enh.pix"),
        dbic=[2, 3, 4,],
        fili_ref=os.path.join(input_dir,"4922_tm_12aug1990_enh.pix"),
        dbic_ref=[2, 3, 4],
        filo=change_image,
        algo="INTENRATIO")

# Extract change as polygons
expolras(fili=change_image,
         dbic=[4],
         thrtype="PixelValue",
         tval=[90],
         areaval=[200],
         filo=change_image)

 

On line 15 filo contains the raster data from after change occurred which is the 1999 dataset. Dbic is set to all six bands within filo. Fili_ref contains the raster data from before the change occurred, which is the 1990 dataset. Dbic_ref is also set to include all 6 bands. On line 19 filo is set as change_image. The output of this algorithm is a file (change_image.pix) with four raster channels. The first channel is the normalized intensity for the new image. The second channel is the normalized intensity for the reference image. The third channel is the absolute value of the intensity ratio (written in decibels) between the new and reference intensity. The fourth channel stores the normalized change, as percentiles. We will be using the fourth channel to create our change polygons.

We can now run the expolras algorithm to extract the change polygons. We will be specifying a threshold value and any pixels above that value will be extracted into the polygons. The input image will be change_image. Use output channel 4, Ranked Changes in Percentile, as the input band. On line 25 and 26 the threshold type is set to PixelValue and then the threshold value is set to 90. Since this tutorial is focusing on deforestation we are interested in large areas of change. We need to specify the minimum area, so we do not extract small clusters of change pixels. Tval is set to 200. We want to save the polygons back to change_image so set it as filo.

Depending on whether you will be doing further analysis with this vector file you can also choose to export it out to a separate file. The vector file includes information on the size of each change area which may be required in other analysis. In our case these vectors will be used to create a binary raster and then an area report.

Polygons showing change from 1990-1999

3. Create binary raster from polygons

Next, we need to create a binary raster from the polygons. This binary raster will be required to run an area report.

 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
# Convert the polygons to a bitmap
poly2bit(fili=change_image,
         dbvs=[2],
         filo=change_image)

# Add a new channel
pcimod(file=change_image,
       pciop="ADD",
       pcival=[1,0,0,0])

# Encode bitmap into an empty image channel
bit_encode(file=change_image,
           dbib=[3],
           valu=[1],
           dboc=[5])

 

The first step for this section will be to convert the vectors to a bitmap. This is done using the poly2bit algorithm. On lines 31 to 33 the input file, input vector segment and output file are set for this algorithm. Again, we will be saving back to change_image

Next, create a new 8-bit channel in input_image. This channel will become the binary raster. On lines 36 to 38 pcimod is used to add this one channel.

In order to now create your binary raster you will need to encode the raster into the empty raster channel you just created. The map algorithm will be used to encode, although you will remember from the import section that this algorithm was imported as bit_encode. For this algorithm you will need to set the input file (file) and input bitmap layer (dbib). The third parameter, valu, specifies that value that will be assigned to the raster pixels below the bitmap. In our case we want this value to be 1. The final parameter is the output image channel. We will set this parameter to the output channel that we created with pcimod, which is 5. Channel 5 will now be a binary raster where the change areas are pixel value 1 and the unchanged areas are pixel value 0.

Binary Raster: Blue sections are change areas

4. Generate area report

The binary raster can now be used in the generation of an area report. There is a report class in the pci.nspio module. You can use this class to save the report information from the areareport to a text file.

 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
# Create report file and generate the area report of changed areas
report_file = os.path.join(working_dir, 'Report.txt')

try:
    Report.clear()
    enableDefaultReport(report_file)
    areareport(file=change_image,
               dbic=[5],
               units="Square Kilometers")
finally:
    enableDefaultReport('term')

 

On lines 47 you will need to set the output path and filename of the report file. Using a python try and except statement you can set up the reporting. On line 50, call the Report.clear() to clear any report information that may still be in memory. You will then need to make sure that the report is being generated in the file instead of the terminal window. Line 51 enables the report implementation and sends the report information to report_file.

With the reporting to file enabled, call the areareport algorithm. The input file will be change_image and we will use the new channel 5 as the input channel. The units for this image will be Square Kilometers. Lastly, add the finally statement and close the report file by calling enableDefaultReport('term').

A report will then be generated which shows the percentage and total areas of change. From the report below you can see that 5.82% of the total area or 824.78 km2 showed significant change from 1990 to 1999.

Have more questions? Submit a request

Comments

Powered by Zendesk