Integrating Landsat, ICESat and ALOS PALSAR for Regional Scale Vegetation Structure Assessment - Source Data

This data, and the following code, allows the reproduction of results from Integrating Landsat, ICESat and ALOS PALSAR for Regional Scale Vegetation Structure Assessment   . All computation was completed using open source software including GDAL, RSGISLib, QGIS and SAGE.

----------injune_2009_fpc_hh_hv_lee--------------

Coregistered FPC and Radar Image over Injune, Australia

Three band geotiff file. Band 1 is foliage projective cover (FPC) downloaded from http://www.derm.qld.gov.au/services_resources/item_details.php?item_id=34767 and bands 2 and 3 are ALOS PALSAR FBD HH and HV data downloaded from:

http://www.eorc.jaxa.jp/ALOS/en/kc_mosaic/kc_mosaic.htm

The original data are provided by JAXA as the ALOS sample product. © JAXA, METI

This is the file used as input to the segmentation below:

----------injune_2009_fpc_hh_hv_clumps_elim_final--------------

Coregistered FPC and Radar Image over Injune, Australia was segmented using the open source  RSGisLib. Input file was the following rsgis xml commands:

rsgis:command algor="imageutils" option="stretch" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_lee.envi" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_stretched.img" ignorezeros="yes" stretch="LinearStdDev" stddev="2" format="HFA"

rsgis:command algor="imagecalc" option="kmeanscentres" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_stretched.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters" numclusters="30" maxiterations="200" degreeofchange="0.25" subsample="10" initmethod="diagonal_range_attach"

rsgis:command algor="segmentation" option="labelsfromclusters" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_stretched.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters.img" clusters="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters.gmtxt" ignorezeros="yes" format="HFA" proj="IMAGE"

rsgis:command algor="segmentation" option="elimsinglepxls" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_stretched.img" clumps="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters.img" temp="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters_singlepxls_tmp.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters_nosinglepxls.img" ignorezeros="yes" format="HFA" proj="IMAGE"

rsgis:command algor="segmentation" option="clump" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_clusters_nosinglepxls.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps.img" nodata="0" format="HFA" inmemory="no" proj="IMAGE"

rsgis:command algor="segmentation" option="rmsmallclumpsstepwise" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_stretched.img" clumps="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim.img" minsize="250" maxspectraldist="200000" format="HFA" inmemory="no" proj="IMAGE"

rsgis:command algor="segmentation" option="relabelclumps" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim_final.img" format="HFA" inmemory="no" proj="IMAGE"

rsgis:command algor="segmentation" option="meanimg" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_lee.envi" clumps="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim_final.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim_mean.img" format="HFA" inmemory="yes" proj="IMAGE"

rsgis:command algor="segmentation" option="randomcolourclumps" image="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim_final.img" output="/injune_project/segmentation/injune_2009_fpc_hh_hv_clumps_elim_clrd.img" format="HFA" inmemory="no" proj="IMAGE"

----------icesat_returns--------------

Subset of GLAS/ICESat GLA14 Release 33 data over the Injune Study Site.

Binary data unpacked and exported as a shapefile

Fields in the dataset correspond to those described in the NSIDC data dictionary except for the SRTM slope field which was calculated from the 1 second SRTM DEM. 

Original data provided by NSIDC.

Zwally, H.J., R. Schutz, C. Bentley, J. Bufton, T. Herring, J. Minster, J. Spinhirne, R. Thomas. 2003, updated current year.GLAS/ICESat L2 Global Land Surface Altimetry Data V018, 15 October to 18 November 2003. Boulder, CO: National Snow and Ice Data Center. Digital media.

----------clustered_segments_with_icesat_heights--------------

Final output shapefile of vegetation heights produced by vectorizing the segmentation image. Segmentation was then clustered using computeClusters.py code below and height attributes were added using the method outlined in the paper, and by the code in the segment_waveform_plots section.

################## computeClusters.py #######################

# Final Clustering Code - reads statistics from computeZonalCovariance

from numpy import load, save

from scipy.spatial.distance import pdist

from scipy.cluster.hierarchy import linkage, fcluster

"""

This script used the statistics from computeZonalCovariance.py

to perform clustering on the segments using the scipy.cluster functions

The output clusterArray maps the input segments to the apprpriate cluster

"""

numClusters=40

# Read in Segment Statistics

print 'Reading in statistics file'

segmentStats=load('segmentStats.npz')

# Compute Compressed Distance Matrix

print 'Computing distance matrix'

D = pdist(segmentStats['meanMatrix'], 'mahalanobis', VI=segmentStats['covMatrix'])

save('distanceMatrix',D)

# Compute Linkage Function

print 'Computing linkage function'

L=linkage(D, method='complete')

save('linkageMatrix',L)

# Compute Flat Clusters

print 'Building flat clusters'

clusterIDX=fcluster(L, numClusters, criterion='maxclust')

save('clusterArray',clusterIDX)

print 'Made ' + str(max(clusterIDX)) + ' clusters'

################## computeZonalCovariance.py #######################

# Little bit of code to compute statistics for mahalanobis distance clustering

from osgeo import gdal

from numpy import unique, mean, cov, sum, savez, array

from scipy import zeros, nonzero

from scipy.linalg import pinv

from joblib import Parallel, delayed

"""

This script computes the zonal covariance based on a segmented raster image

and the original values raster. The output file is used to cluster the segments

using computeClusters.py

"""

clumpedRaster='injune_2009_fpc_hh_hv_1000_clumps_elim_final.img'

valuesRaster='injune_2009_fpc_hh_hv_lee.tif'

# Open Categories Raster

dsC = gdal.Open(clumpedRaster)

categoryRaster = dsC.GetRasterBand(1).ReadAsArray().flatten(1)

categoryBins=unique(categoryRaster.flatten(1))

# Open Values Raster

dsV = gdal.Open(valuesRaster)

numBands=dsV.RasterCount

numValues=categoryRaster.size

numCategories=categoryBins.size

# Read in all Bands

valuesRaster=zeros((numBands,numValues))

for band in range(numBands):

valuesRaster[band] = dsV.GetRasterBand(band+1).ReadAsArray().flatten(1)

# Compute Statistics

def computeStats(categoryID):

print numCategories-categoryID

pixelIDX=(categoryRaster==categoryBins[categoryID])

if sum(pixelIDX)>3: # Make sure we can calculate some stats

return [categoryBins[categoryID],valuesRaster[:,pixelIDX].mean(axis=1),pinv(cov(valuesRaster[:,pixelIDX],rowvar = True))]

else:

return [0, array([0,0,0]), array([[0,0,0],[0,0,0],[0,0,0]])]

allStats=Parallel(n_jobs=8)(delayed(computeStats)(categoryID) for categoryID in range(numCategories))

categoryMatrix, meanMatrix, covMatrix = zip(*allStats)

# Remove Zeros

goodDataIDX=(array(categoryMatrix) > 0)

categoryMatrix=array(categoryMatrix)[goodDataIDX]

meanMatrix=array(meanMatrix)[goodDataIDX]

covMatrix=array(covMatrix)[goodDataIDX]

# Write data to text file

savez('segmentStats',goodDataIDX=goodDataIDX,categoryMatrix=categoryMatrix,meanMatrix=meanMatrix,covMatrix=covMatrix)

----------segment_waveform_plots--------------

Vegetation vertical profile plots representing the mean ICESat waveforms captured within each cluster of similar segments. CSV file has summary statistics. Figures on plots are, from top, the cluster ID, Height to the greatest vegetation density, height where 50% of the returns have been intercepted, height where 95% of the returns have been intercepted, number of ICESat returns within the cluster, estimated cover in percent and ground surface standard deviation in meters.

Figures were generated using the following python code:

# Code to generate the aggregated pulses in the clustered segments

"""

Note: THIS IS A SAGE SCRIPT so you'll need to either import sage, or run it within a sage notebook.

If you want to run this you will probably need some more information than my sketchy 

documentation of this messy code so drop me an email at [e-mail]

Reads a text file containing all the icesat waveforms within each segment

Need to preprocess the spatially joined dbf using:

dbview --browse --delimiter=, gla14_r28_qld_23sProject_Spa.dbf | awk 'BEGIN { FS = "," };{ gsub(" ",""); print $33,$35,$36,$37,$38,$40,$41,$42,$44,$45,$46,$48,$49,$50,$52,$53,$54,$56,$57,$58,$60 >> $62 ".txt" }'

"""

import os

import csv

import sys

from scipy import arange, argmin, argmax, cumsum, loadtxt, nonzero, zeros

from scipy.optimize import fmin

from scipy.signal import deconvolve

from scipy.linalg import pinv

from numpy import dot,finfo

# Get the filenames

reDataDir='/home/pete/share/Dropbox/injune_ICESat/1000Segment/intersect/'

plotDataDir='/home/pete/share/Dropbox/injune_ICESat/1000Segment/plots/'

filenames=os.listdir(reDataDir)

eps = finfo(float).eps

def ellipseArea(majorAxis,eccentricity):

    return 3.14159*majorAxis*sqrt(majorAxis^2-eccentricity^2*majorAxis^2)

def gaussian(x,offset,amplitude,sigma):

    return amplitude*exp(-(x-offset)^2/(2*(sigma/2.99792458+eps)^2)) # NOTE - 2.99792458 is conversion for release 33 format reverting to nanoseconds

def iceSatWaveform(x, gaussianFits):

    # Find the "ground pulse"

    minWaveform=argmin(gaussianFits[[0,3,6,9,12,15]])

    minOffset=gaussianFits[minWaveform]

    # Recover the ground pulse area

    groundPulse=gaussian(x,gaussianFits[minWaveform]-minOffset,\

    gaussianFits[minWaveform+1],gaussianFits[minWaveform+2])

    # Set the amplitude of the ground pulse to zero before building the vegetation waveform

    gaussianFits[minWaveform+1]=0

    waveform= \

    gaussian(x,gaussianFits[0]-minOffset,gaussianFits[1],gaussianFits[2])+\

    gaussian(x,gaussianFits[3]-minOffset,gaussianFits[4],gaussianFits[5])+\

    gaussian(x,gaussianFits[6]-minOffset,gaussianFits[7],gaussianFits[8])+\

    gaussian(x,gaussianFits[9]-minOffset,gaussianFits[10],gaussianFits[11])+\

    gaussian(x,gaussianFits[12]-minOffset,gaussianFits[13],gaussianFits[14])+\

    gaussian(x,gaussianFits[15]-minOffset,gaussianFits[16],gaussianFits[17])

    return [waveform,groundPulse]

def normalisedIceSatWaveform(x, pulseParameters):

    # This normalises the returned waveform by intensity

    rawWaveforms=iceSatWaveform(x,pulseParameters[3:21])

    return rawWaveforms/pulseParameters[0] 

    #*ellipseArea(pulseParameters[2],pulseParameters[1])

def slopeCorrectedIceSatWaveform(x,uncorrectedWaveform,groundSlope,pulseDiameter):

    # This function generates an approximation of the effect of the slope on the retrieved waveform

    zeroIDX=argmax(1 * (x==0))

    signalAtGround=uncorrectedWaveform[zeroIDX]

    fwhmSlopeSignal=pulseDiameter * 100 * 3* sin(groundSlope/180*3.14159+0.000001) #  i_tpmajoraxis_avg is in meters

    slopeSignal=signalAtGround*exp(-(x^2)/(2*(fwhmSlopeSignal)^2))

    correctedWaveform=uncorrectedWaveform-slopeSignal

    correctedWaveform[correctedWaveform<0]=0

    return correctedWaveform

# Generate the height bins to reconstruct the waveform

heightBins=(arange(7001)-2000.0)/100

# Open the summary csv file

reSummaryFile=open(plotDataDir+"summaryRE.csv",'w')

reSummaryOutput=csv.writer(reSummaryFile)

# Loop through all the filesaggregatedWaveformArray=

for reNum in range(len(filenames)):

#for reNum in range(2):

    # Import the satellite data

    satData=loadtxt(reDataDir+filenames[reNum],delimiter=' ')

    # Check if there is more than one VALID pulse on a ground slope of less than 5 degrees

    if (satData.ndim > 1) and not((satData[:,21]>5).all() or (satData[:,0]10000).all()):

        # Filter out crap points

        satData=satData[(satData[:,0]<10000).nonzero()]

        satData=satData[(satData[:,0]>1).nonzero()]

        satData=satData[(satData[:,21]<5).nonzero()]

        # Initialise arrays to hold the cumulative pulses

        numHits=satData.shape[0]

        aggregatedWaveformArray=zeros((numHits,heightBins.size))

        aggregatedGroundArray=zeros((numHits,heightBins.size))

        # Initalise Arrays to compute the relative reflectances

        reflectanceTotal=satData[:,22]

        areaGroundVeg=zeros([satData[:,22].size,2])

        # Kludge for if there is no foliage waveform

        #aggregatedWaveformArray[:,0]=0.0000001

        #aggregatedGroundArray[:,0]=0.0001

        # Loop through all the pulses in the RE

        for i in range(numHits):

            splitWaveforms=normalisedIceSatWaveform(heightBins,satData[i])

            slopeCorrectedVegetationWaveform=slopeCorrectedIceSatWaveform(heightBins,splitWaveforms[0],satData[i,21],satData[i,2])

            aggregatedWaveformArray[i]=slopeCorrectedVegetationWaveform #splitWaveforms[0]

            aggregatedGroundArray[i]=splitWaveforms[1]

            areaGroundVeg[i]=[sum(splitWaveforms[0]),sum(splitWaveforms[1])]/sum(splitWaveforms)

       # Reflectance Ratio Calculation

        reflectances=dot(pinv(areaGroundVeg),reflectanceTotal)

        print reflectances

        # Compute the foliage profile

        aggregatedWaveform=aggregatedWaveformArray.mean(axis=0)

        aggregatedWaveformStdev=aggregatedWaveformArray.std(axis=0)/10

        aggregatedGround=aggregatedGroundArray.mean(axis=0)

        #foliageProfile=cumsum(aggregatedWaveform)/(sum(aggregatedGround+aggregatedWaveform))

        foliageProfile=reflectances[0]*cumsum(aggregatedWaveform)/(reflectances[1]*\

        sum(aggregatedGround)+reflectances[0]*sum(aggregatedWaveform))

        scaledAggregatedWaveform=aggregatedWaveform/aggregatedWaveform.max()*foliageProfile.max()

        scaledAggregatedWaveformStdevPlus=(aggregatedWaveform+aggregatedWaveformStdev)/aggregatedWaveform.max()*foliageProfile.max()

        scaledAggregatedWaveformStdevMinus=(aggregatedWaveform-aggregatedWaveformStdev)/aggregatedWaveform.max()*foliageProfile.max()

        # Some structure metrics

        fpc=foliageProfile.max()

        heightMode=0

        heightAvg=0

        height95=0

        if fpc>0.01:

            heightMode=heightBins[argmax(aggregatedWaveform)]

            heightAvg=heightBins[nonzero(foliageProfile > fpc*0.5)[0][0]]

            height95=heightBins[nonzero(foliageProfile > fpc*0.95)[0][0]]

        # Recover the sigma of the ground pulse as an estimate of suface roughness/groundcover

        def gaussianFit(gaussianParams):

            return sum((aggregatedGround-\

            gaussian(heightBins,gaussianParams[0],gaussianParams[1],gaussianParams[2]))**2)

        groundFitParameters = fmin(gaussianFit, [0.0,aggregatedGround.max(),1])

        # Write structural data to the summary file

        reName=os.path.splitext(filenames[reNum])[0]

        reSummaryOutput.writerow([reName,heightMode,heightAvg,height95,numHits,round(fpc,2),round(groundFitParameters[2],2)])

        # Waveform Plot

        foliageProfilePlot=list_plot(zip(scaledAggregatedWaveform[2000:5201],\

        heightBins[2000:5201]),rgbcolor=(0.3,1.0,0.3),gridlines='automatic',frame=true,figsize=[4,6],\

        plotjoined=True,thickness=3,xmin=0,xmax=fpc,ymin=0,ymax=30)

        foliageProfilePlot=foliageProfilePlot+list_plot(zip(scaledAggregatedWaveformStdevPlus[2000:5201],\

        heightBins[2000:5201]),rgbcolor=(1.0,0.6,0.0),thickness=2,linestyle=":",plotjoined=True)

        foliageProfilePlot=foliageProfilePlot+list_plot(zip(scaledAggregatedWaveformStdevMinus[2000:5201],\

        heightBins[2000:5201]),rgbcolor=(1.0,0.6,0.0),thickness=2,linestyle=":",plotjoined=True)

        foliageProfilePlot=foliageProfilePlot+list_plot(zip(foliageProfile[2000:5201],heightBins[2000:5201]),\

        rgbcolor=(0.4,0.5,0.0),plotjoined=True,thickness=3)

        foliageProfilePlot=foliageProfilePlot+text(reName, (fpc/2,29), fontsize=14, color='darkred')+\

            text(heightMode, (fpc/2,27))+\

            text(heightAvg, (fpc/2,26))+\

            text(height95, (fpc/2,25))+\

            text(numHits, (fpc/2,23))+\

            text(round(fpc*100,1), (fpc/2,22))+\

            text(round(2*groundFitParameters[2]/2.99792458,1), (fpc/2,21))

        foliageProfilePlot.axes_labels(['Cover %','Height (m)']) 

        # Export the foliage profile plot

        foliageProfilePlot.save(plotDataDir+reName+".png",gridlines='automatic',frame=true,figsize=[4,6])

# Close the summary file

reSummaryFile.close()