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

• Main Author: Scarth, Peter
• Format: Dataset
• Language: unknown
• Published: figshare 2012
• Subjects: Physical Geography; Environmental Science; Ecology; FOS Biological sciences; National Snow and Ice Data Center
• Online Access: https://dx.doi.org/10.6084/m9.figshare.94245.v1; https://figshare.com/articles/dataset/Integrating_Landsat,_ICESat_and_ALOS_PALSAR_for_Regional_Scale_Vegetation_Structure_Assessment_-_Source_Data/94245/1

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()