Examples of pkfilter

Filtering in spatial domain

Filter input.tif with morphological dilation filter. Use a circular kernel (instead of rectangular) of size 3x3.

pkfilter -i input.tif -o filter.tif -dx 3 -dy 3 -f dilate -circ

Similar to previous example, but consider only values of 255 for filtering operation. Typical usage: dilate cloud values in input image that are flagged as 255

pkfilter -i input.tif -o filter.tif -dx 3 -dy 3 -class 255 -f dilate -circ

Filtering in spectral/temporal domain

Applications with moving window

Calculate the median value for each pixel, calculated on a moving window of width 3 (-dz 3) over all input bands. The output raster dataset will contain as many bands as the input raster dataset.

pkfilter -i input.tif -o filter_stdev.tif -dz 3 -f median

Applications of statistical functions in spectral domain

Calculate the standard deviation for each pixel, calculated on all input bands. The output raster dataset will contain a single band only, no moving window is used (-dz 1).

pkfilter -i input.tif -o filter_stdev.tif -dz 1 -f stdev

"Smooth" (interpolate) nodata in spectral/temporal domain (-dz 1), using a linear interpolation. The following interpolation types are supported: akima (default), linear, polynomial, cspline, cspline_periodic, akima_periodic (please check gsl page for more information on the interpolation types).

pkfilter -i input.tif -o smoothed.tif -dz 1 -f smoothnodata -interp linear

Filter with spectral response functions

The following two examples show how to use pkfilter for spectral filtering a high dimensional input (hyperspectral image) to a lower dimensional output (multi-spectral image). Notice that the input wavelenghts must be provided as -win value1 -win value2 -win value3 ... To save typing, we assume the input wavelengths are listed in a text file wavelengths.txt (single column ASCII file with all wavelenghts listed in nanometer).

cat wavelengths.txt | while read W;do echo " -win $W";done

In the first example, the hyperspectral image is filtered with a spectral response function. For each spectral response function provided, a separate output band is created. The spectral response function(s) must be listed in two column ASCII file(s) with the wavelengths and response listed in the first and second column respectively. The response functions can but must not be normalized (this is taken care of by the filter utility).

In this example, the input is a hyperspectral image with N>>1 spectral wavelengths (bands). The output is a multispectral image with 3 bands, where a spectral response function is provided for each output band.

pkfilter -i hyperspectral.tif -o multispectral.tif -srf srf1.txt -srf srf2.txt -srf srf3.txt $(cat wavelengths.txt | while read W;do echo " -win $W";done)

The next example is similar to the previous. Instead of providing a spectral response function for each output band, you can also provide the center wavelengths and full width half max values. Here, a three band (red, green, blue) output image is produced.

pkfilter -i hyperspectral.tif -o multispectral.tif -wout 650 -wout 510 -wout 475 -fwhm 50 -fwhm 50 -fwhm 50 $(cat wavelengths.txt | while read W;do echo " -win $W";done)

Applying the Savitzky-Golay filter to reconstruct a time-series data set

The following example reconstructs a time-series data set based on the Savitzky–Golay filter, as suggested by J. Chen 2004.

Input is a multi-band (Byte) raster dataset containing a noisy time series, e.g., a normalized difference vegetation index (NDVI) with cloud contaminated (low) NDVI values. Output is the reconstructed time-series data set, which approaches the upper NDVI enveloppe. Please refer to J. Chen 2004 for more details.

Preparation

Create a long-term change trend fitting (lta.tif) using the Savitzky-Golay filter. Choose the number of lefward (past) and rightward (future) data points (e.g., 7). The order of smoothing polynomial in the Savitzky-Golay filter is set to 2.

pkfilter -i input.tif -o lta.tif -f savgolay -nl 7 -nr 7 -m 2 -pad replicate

Initialization

Decrease the number of lefward (past) and rightward (future) data points to 4. The order of smoothing polynomial in the Savitzky-Golay filter is set to 6. Valid data are from 0 to 250 (nodata is set to 250).

pkcomposite -i lta.tif -i input.tif -o iter1.tif -cr maxallbands -max 250 -dstnodata 250

pkfilter -i iter1.tif -o savgolay1.tif -f savgolay -nl 4 -nr 4 -m 6 -pad replicate

Iterative process (repeat, e.g., 5-10 times)

pkcomposite -i input.tif -i savgolay1.tif -o iter2.tif -cr maxallbands -max 250 -dstnodata

pkfilter -i iter2.tif -o savgolay2.tif -f savgolay -nl 4 -nr 4 -m 6 -pad replicate

pkcomposite -i input.tif -i savgolay2.tif -o iter3.tif -cr maxallbands -max 250 -dstnodata

pkfilter -i iter3.tif -o savgolay3.tif -f savgolay -nl 4 -nr 4 -m 6 -pad replicate

etc.

Output savgolay<n>.tif is the reconstructed time-series data set.