10.6084/m9.figshare.94245.v1 Peter Scarth Peter Scarth Integrating Landsat, ICESat and ALOS PALSAR for Regional Scale Vegetation Structure Assessment - Source Data figshare 2012 ALOS ICESat Landsat Lidar OBIA segmentation Vegetation Structure Physical Geography Environmental Science Ecology 2012-08-11 12:34:25 Dataset https://figshare.com/articles/dataset/Integrating_Landsat_ICESat_and_ALOS_PALSAR_for_Regional_Scale_Vegetation_Structure_Assessment_-_Source_Data/94245 <p> </p> <p><a></a><a></a> This data, and the following code, allows the reproduction of results from <a href="http://dx.doi.org/10.6084/m9.figshare.94112">Integrating Landsat, ICESat and ALOS PALSAR for Regional Scale Vegetation Structure Assessment </a>  . All computation was completed using open source software including <a href="http://www.gdal.org/">GDAL</a>, <a href="http://www.rsgislib.org/">RSGISLib</a>, <a href="http://www.qgis.org/">QGIS</a> and <a href="http://www.sagemath.org/">SAGE</a>.</p> <p><br><br></p> <p><strong>----------injune_2009_fpc_hh_hv_lee--------------</strong></p> <p>Coregistered FPC and Radar Image over Injune, Australia</p> <p>Three band geotiff file. Band 1 is foliage projective cover (FPC) downloaded from <a href="http://www.derm.qld.gov.au/services_resources/item_details.php?item_id=34767">http://www.derm.qld.gov.au/services_resources/item_details.php?item_id=34767</a> and bands 2 and 3 are ALOS PALSAR FBD HH and HV data downloaded from:</p> <p><a href="http://www.eorc.jaxa.jp/ALOS/en/kc_mosaic/kc_mosaic.htm">http://www.eorc.jaxa.jp/ALOS/en/kc_mosaic/kc_mosaic.htm</a></p> <p>The original data are provided by JAXA as the ALOS sample product. © JAXA, METI</p> <p>This is the file used as input to the segmentation below:</p> <p><strong>----------injune_2009_fpc_hh_hv_clumps_elim_final--------------</strong></p> <p>Coregistered FPC and Radar Image over Injune, Australia was segmented using the open source  <a href="http://www.rsgislib.org/">RSGisLib</a>. Input file was the following rsgis xml commands:</p> <p>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" </p> <p>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"</p> <p> 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"</p> <p>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"</p> <p>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" </p> <p>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"</p> <p>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"</p> <p>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" </p> <p>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"</p> <p><br><br></p> <p><strong>----------icesat_returns--------------</strong></p> <p>Subset of <a href="http://nsidc.org/data/gla14.html">GLAS/ICESat GLA14</a> Release 33 data over the Injune Study Site.</p> <p>Binary data unpacked and exported as a shapefile</p> <p>Fields in the dataset correspond to those described in the<a href="http://nsidc.org/data/docs/daac/glas_altimetry/gla14_records.html"> NSIDC data dictionary</a> except for the SRTM slope field which was calculated from the 1 second SRTM DEM. </p> <p>Original data provided by NSIDC.</p> <p><em>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.</em></p> <p><br><br></p> <p><strong>----------clustered_segments_with_icesat_heights--------------</strong></p> <p><a></a><a></a> 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.</p> <p> </p> <p><a></a> <strong>################## computeClusters.py #######################</strong></p> <p> </p> <p># Final Clustering Code - reads statistics from computeZonalCovariance </p> <p> </p> <p>from numpy import load, save </p> <p>from scipy.spatial.distance import pdist </p> <p>from scipy.cluster.hierarchy import linkage, fcluster </p> <p> </p> <p>""" </p> <p>This script used the statistics from computeZonalCovariance.py </p> <p>to perform clustering on the segments using the scipy.cluster functions </p> <p>The output clusterArray maps the input segments to the apprpriate cluster </p> <p>""" </p> <p> </p> <p>numClusters=40 </p> <p> </p> <p># Read in Segment Statistics </p> <p>print 'Reading in statistics file' </p> <p>segmentStats=load('segmentStats.npz') </p> <p> </p> <p># Compute Compressed Distance Matrix </p> <p>print 'Computing distance matrix' </p> <p>D = pdist(segmentStats['meanMatrix'], 'mahalanobis', VI=segmentStats['covMatrix']) </p> <p>save('distanceMatrix',D) </p> <p> </p> <p># Compute Linkage Function </p> <p>print 'Computing linkage function' </p> <p>L=linkage(D, method='complete') </p> <p>save('linkageMatrix',L) </p> <p> </p> <p># Compute Flat Clusters </p> <p>print 'Building flat clusters' </p> <p>clusterIDX=fcluster(L, numClusters, criterion='maxclust') </p> <p>save('clusterArray',clusterIDX) </p> <p> </p> <p>print 'Made ' + str(max(clusterIDX)) + ' clusters'</p> <p><br><br></p> <p><strong>################## computeZonalCovariance.py #######################</strong></p> <p> </p> <p># Little bit of code to compute statistics for mahalanobis distance clustering </p> <p> </p> <p>from osgeo import gdal </p> <p>from numpy import unique, mean, cov, sum, savez, array </p> <p>from scipy import zeros, nonzero </p> <p>from scipy.linalg import pinv </p> <p>from joblib import Parallel, delayed </p> <p> </p> <p>""" </p> <p>This script computes the zonal covariance based on a segmented raster image </p> <p>and the original values raster. The output file is used to cluster the segments </p> <p>using computeClusters.py </p> <p>""" </p> <p> </p> <p>clumpedRaster='injune_2009_fpc_hh_hv_1000_clumps_elim_final.img' </p> <p>valuesRaster='injune_2009_fpc_hh_hv_lee.tif' </p> <p> </p> <p># Open Categories Raster </p> <p>dsC = gdal.Open(clumpedRaster) </p> <p>categoryRaster = dsC.GetRasterBand(1).ReadAsArray().flatten(1) </p> <p>categoryBins=unique(categoryRaster.flatten(1)) </p> <p> </p> <p> </p> <p># Open Values Raster </p> <p>dsV = gdal.Open(valuesRaster) </p> <p>numBands=dsV.RasterCount </p> <p>numValues=categoryRaster.size </p> <p>numCategories=categoryBins.size </p> <p> </p> <p># Read in all Bands </p> <p>valuesRaster=zeros((numBands,numValues)) </p> <p> </p> <p>for band in range(numBands): </p> <p> valuesRaster[band] = dsV.GetRasterBand(band+1).ReadAsArray().flatten(1) </p> <p> </p> <p># Compute Statistics </p> <p> </p> <p>def computeStats(categoryID): </p> <p> print numCategories-categoryID </p> <p> pixelIDX=(categoryRaster==categoryBins[categoryID]) </p> <p> if sum(pixelIDX)>3: # Make sure we can calculate some stats </p> <p> return [categoryBins[categoryID],valuesRaster[:,pixelIDX].mean(axis=1),pinv(cov(valuesRaster[:,pixelIDX],rowvar = True))] </p> <p> else: </p> <p> return [0, array([0,0,0]), array([[0,0,0],[0,0,0],[0,0,0]])] </p> <p> </p> <p>allStats=Parallel(n_jobs=8)(delayed(computeStats)(categoryID) for categoryID in range(numCategories)) </p> <p> </p> <p>categoryMatrix, meanMatrix, covMatrix = zip(*allStats) </p> <p> </p> <p># Remove Zeros </p> <p>goodDataIDX=(array(categoryMatrix) > 0) </p> <p>categoryMatrix=array(categoryMatrix)[goodDataIDX] </p> <p>meanMatrix=array(meanMatrix)[goodDataIDX] </p> <p>covMatrix=array(covMatrix)[goodDataIDX] </p> <p> </p> <p> </p> <p># Write data to text file </p> <p>savez('segmentStats',goodDataIDX=goodDataIDX,categoryMatrix=categoryMatrix,meanMatrix=meanMatrix,covMatrix=covMatrix) </p> <p> </p> <p><br><br></p> <p><strong>----------segment_waveform_plots--------------</strong></p> <p>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.</p> <p>Figures were generated using the following python code:</p> <p><br><br></p> <p># Code to generate the aggregated pulses in the clustered segments</p> <p> </p> <p>"""</p> <p>Note: THIS IS A SAGE SCRIPT so you'll need to either import sage, or run it within a sage notebook.</p> <p>If you want to run this you will probably need some more information than my sketchy </p> <p>documentation of this messy code so drop me an email at p.scarth@uq.edu.au</p> <p> </p> <p>Reads a text file containing all the icesat waveforms within each segment</p> <p>Need to preprocess the spatially joined dbf using:</p> <p>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" }'</p> <p>"""</p> <p>import os</p> <p>import csv</p> <p>import sys</p> <p>from scipy import arange, argmin, argmax, cumsum, loadtxt, nonzero, zeros</p> <p>from scipy.optimize import fmin</p> <p>from scipy.signal import deconvolve</p> <p>from scipy.linalg import pinv</p> <p>from numpy import dot,finfo</p> <p> </p> <p># Get the filenames</p> <p>reDataDir='/home/pete/share/Dropbox/injune_ICESat/1000Segment/intersect/'</p> <p>plotDataDir='/home/pete/share/Dropbox/injune_ICESat/1000Segment/plots/'</p> <p>filenames=os.listdir(reDataDir)</p> <p> </p> <p>eps = finfo(float).eps</p> <p> </p> <p>def ellipseArea(majorAxis,eccentricity):</p> <p>    return 3.14159*majorAxis*sqrt(majorAxis^2-eccentricity^2*majorAxis^2)</p> <p> </p> <p>def gaussian(x,offset,amplitude,sigma):</p> <p>    return amplitude*exp(-(x-offset)^2/(2*(sigma/2.99792458+eps)^2)) # NOTE - 2.99792458 is conversion for release 33 format reverting to nanoseconds</p> <p> </p> <p>def iceSatWaveform(x, gaussianFits):</p> <p>    # Find the "ground pulse"</p> <p>    minWaveform=argmin(gaussianFits[[0,3,6,9,12,15]])</p> <p>    minOffset=gaussianFits[minWaveform]</p> <p>    </p> <p>    # Recover the ground pulse area</p> <p>    groundPulse=gaussian(x,gaussianFits[minWaveform]-minOffset,\</p> <p>    gaussianFits[minWaveform+1],gaussianFits[minWaveform+2])</p> <p>    </p> <p>    # Set the amplitude of the ground pulse to zero before building the vegetation waveform</p> <p>    gaussianFits[minWaveform+1]=0</p> <p>    waveform= \</p> <p>    gaussian(x,gaussianFits[0]-minOffset,gaussianFits[1],gaussianFits[2])+\</p> <p>    gaussian(x,gaussianFits[3]-minOffset,gaussianFits[4],gaussianFits[5])+\</p> <p>    gaussian(x,gaussianFits[6]-minOffset,gaussianFits[7],gaussianFits[8])+\</p> <p>    gaussian(x,gaussianFits[9]-minOffset,gaussianFits[10],gaussianFits[11])+\</p> <p>    gaussian(x,gaussianFits[12]-minOffset,gaussianFits[13],gaussianFits[14])+\</p> <p>    gaussian(x,gaussianFits[15]-minOffset,gaussianFits[16],gaussianFits[17])</p> <p>            </p> <p>    return [waveform,groundPulse]</p> <p> </p> <p> </p> <p>def normalisedIceSatWaveform(x, pulseParameters):</p> <p>    # This normalises the returned waveform by intensity</p> <p>    rawWaveforms=iceSatWaveform(x,pulseParameters[3:21])</p> <p>    return rawWaveforms/pulseParameters[0] </p> <p>    #*ellipseArea(pulseParameters[2],pulseParameters[1])</p> <p>    </p> <p> </p> <p>def slopeCorrectedIceSatWaveform(x,uncorrectedWaveform,groundSlope,pulseDiameter):</p> <p>    # This function generates an approximation of the effect of the slope on the retrieved waveform</p> <p>    zeroIDX=argmax(1 * (x==0))</p> <p>    signalAtGround=uncorrectedWaveform[zeroIDX]</p> <p>    fwhmSlopeSignal=pulseDiameter * 100 * 3* sin(groundSlope/180*3.14159+0.000001) #  i_tpmajoraxis_avg is in meters</p> <p>    slopeSignal=signalAtGround*exp(-(x^2)/(2*(fwhmSlopeSignal)^2))</p> <p>    correctedWaveform=uncorrectedWaveform-slopeSignal</p> <p>    correctedWaveform[correctedWaveform<0]=0</p> <p>    return correctedWaveform</p> <p>    </p> <p> </p> <p># Generate the height bins to reconstruct the waveform</p> <p>heightBins=(arange(7001)-2000.0)/100</p> <p> </p> <p># Open the summary csv file</p> <p>reSummaryFile=open(plotDataDir+"summaryRE.csv",'w')</p> <p>reSummaryOutput=csv.writer(reSummaryFile)</p> <p> </p> <p> </p> <p># Loop through all the filesaggregatedWaveformArray=</p> <p> </p> <p>for reNum in range(len(filenames)):</p> <p>#for reNum in range(2):</p> <p> </p> <p>    # Import the satellite data</p> <p>    </p> <p>    satData=loadtxt(reDataDir+filenames[reNum],delimiter=' ')</p> <p>    </p> <p>    # Check if there is more than one VALID pulse on a ground slope of less than 5 degrees</p> <p>    if (satData.ndim > 1) and not((satData[:,21]>5).all() or (satData[:,0]10000).all()):</p> <p>    </p> <p>        # Filter out crap points</p> <p>        satData=satData[(satData[:,0]<10000).nonzero()]</p> <p>        satData=satData[(satData[:,0]>1).nonzero()]</p> <p>        satData=satData[(satData[:,21]<5).nonzero()]</p> <p>        </p> <p>        # Initialise arrays to hold the cumulative pulses</p> <p>        numHits=satData.shape[0]</p> <p>        aggregatedWaveformArray=zeros((numHits,heightBins.size))</p> <p>        aggregatedGroundArray=zeros((numHits,heightBins.size))</p> <p>       </p> <p>        # Initalise Arrays to compute the relative reflectances</p> <p>        reflectanceTotal=satData[:,22]</p> <p>        areaGroundVeg=zeros([satData[:,22].size,2])</p> <p> </p> <p>       </p> <p>        # Kludge for if there is no foliage waveform</p> <p>        #aggregatedWaveformArray[:,0]=0.0000001</p> <p>        #aggregatedGroundArray[:,0]=0.0001</p> <p>        </p> <p>        # Loop through all the pulses in the RE</p> <p>        for i in range(numHits):</p> <p>            splitWaveforms=normalisedIceSatWaveform(heightBins,satData[i])</p> <p>            slopeCorrectedVegetationWaveform=slopeCorrectedIceSatWaveform(heightBins,splitWaveforms[0],satData[i,21],satData[i,2])</p> <p>            aggregatedWaveformArray[i]=slopeCorrectedVegetationWaveform #splitWaveforms[0]</p> <p>            aggregatedGroundArray[i]=splitWaveforms[1]</p> <p>            areaGroundVeg[i]=[sum(splitWaveforms[0]),sum(splitWaveforms[1])]/sum(splitWaveforms)</p> <p>        </p> <p>     </p> <p>       # Reflectance Ratio Calculation</p> <p>        reflectances=dot(pinv(areaGroundVeg),reflectanceTotal)</p> <p>        print reflectances</p> <p> </p> <p>        # Compute the foliage profile</p> <p>        aggregatedWaveform=aggregatedWaveformArray.mean(axis=0)</p> <p>        aggregatedWaveformStdev=aggregatedWaveformArray.std(axis=0)/10</p> <p>        </p> <p>        aggregatedGround=aggregatedGroundArray.mean(axis=0)</p> <p> </p> <p>        #foliageProfile=cumsum(aggregatedWaveform)/(sum(aggregatedGround+aggregatedWaveform))</p> <p>        foliageProfile=reflectances[0]*cumsum(aggregatedWaveform)/(reflectances[1]*\</p> <p>        sum(aggregatedGround)+reflectances[0]*sum(aggregatedWaveform))</p> <p>        scaledAggregatedWaveform=aggregatedWaveform/aggregatedWaveform.max()*foliageProfile.max()</p> <p>        scaledAggregatedWaveformStdevPlus=(aggregatedWaveform+aggregatedWaveformStdev)/aggregatedWaveform.max()*foliageProfile.max()</p> <p>        scaledAggregatedWaveformStdevMinus=(aggregatedWaveform-aggregatedWaveformStdev)/aggregatedWaveform.max()*foliageProfile.max()</p> <p>        </p> <p>        </p> <p>        # Some structure metrics</p> <p>        fpc=foliageProfile.max()</p> <p>        heightMode=0</p> <p>        heightAvg=0</p> <p>        height95=0</p> <p>        if fpc>0.01:</p> <p>            heightMode=heightBins[argmax(aggregatedWaveform)]</p> <p>            heightAvg=heightBins[nonzero(foliageProfile > fpc*0.5)[0][0]]</p> <p>            height95=heightBins[nonzero(foliageProfile > fpc*0.95)[0][0]]</p> <p>            </p> <p>        # Recover the sigma of the ground pulse as an estimate of suface roughness/groundcover</p> <p>        def gaussianFit(gaussianParams):</p> <p>            return sum((aggregatedGround-\</p> <p>            gaussian(heightBins,gaussianParams[0],gaussianParams[1],gaussianParams[2]))**2)</p> <p>        </p> <p> </p> <p>            </p> <p>        groundFitParameters = fmin(gaussianFit, [0.0,aggregatedGround.max(),1])</p> <p>        </p> <p>        # Write structural data to the summary file</p> <p>        reName=os.path.splitext(filenames[reNum])[0]</p> <p>        reSummaryOutput.writerow([reName,heightMode,heightAvg,height95,numHits,round(fpc,2),round(groundFitParameters[2],2)])</p> <p>     </p> <p>        </p> <p>        </p> <p>        # Waveform Plot</p> <p>        </p> <p>        foliageProfilePlot=list_plot(zip(scaledAggregatedWaveform[2000:5201],\</p> <p>        heightBins[2000:5201]),rgbcolor=(0.3,1.0,0.3),gridlines='automatic',frame=true,figsize=[4,6],\</p> <p>        plotjoined=True,thickness=3,xmin=0,xmax=fpc,ymin=0,ymax=30)</p> <p>        foliageProfilePlot=foliageProfilePlot+list_plot(zip(scaledAggregatedWaveformStdevPlus[2000:5201],\</p> <p>        heightBins[2000:5201]),rgbcolor=(1.0,0.6,0.0),thickness=2,linestyle=":",plotjoined=True)</p> <p>        foliageProfilePlot=foliageProfilePlot+list_plot(zip(scaledAggregatedWaveformStdevMinus[2000:5201],\</p> <p>        heightBins[2000:5201]),rgbcolor=(1.0,0.6,0.0),thickness=2,linestyle=":",plotjoined=True)</p> <p>        foliageProfilePlot=foliageProfilePlot+list_plot(zip(foliageProfile[2000:5201],heightBins[2000:5201]),\</p> <p>        rgbcolor=(0.4,0.5,0.0),plotjoined=True,thickness=3)</p> <p>        foliageProfilePlot=foliageProfilePlot+text(reName, (fpc/2,29), fontsize=14, color='darkred')+\</p> <p>            text(heightMode, (fpc/2,27))+\</p> <p>            text(heightAvg, (fpc/2,26))+\</p> <p>            text(height95, (fpc/2,25))+\</p> <p>            text(numHits, (fpc/2,23))+\</p> <p>            text(round(fpc*100,1), (fpc/2,22))+\</p> <p>            text(round(2*groundFitParameters[2]/2.99792458,1), (fpc/2,21))</p> <p>                </p> <p>        foliageProfilePlot.axes_labels(['Cover %','Height (m)']) </p> <p>      </p> <p> </p> <p>        </p> <p>        # Export the foliage profile plot</p> <p>        foliageProfilePlot.save(plotDataDir+reName+".png",gridlines='automatic',frame=true,figsize=[4,6])</p> <p>        </p> <p># Close the summary file</p> <p>reSummaryFile.close()</p> <p> </p> <p> </p> <p><a></a><a></a><a></a> </p>