regarding DQMeasure : Python analysis tools seem to miss-read unit factor (feet/metre) (possibly a simple python question)

[Sam Hackett]

Hi Guys, wondering if any of you are familiar with DQMeasure and the python analysis tools designed provide some stats from the output (DQMOutputAnalysisHorizontal.py with DQMAnalysisFun_Horizontal.py) (running in Python 2.7.14)

The scripts run for me, however it appears DQMAnalysisFun_Horizontal.py is confusing the unit factor input value (1 = metre, 0.3048 = feet).  I am thinking this is the case because the result .bat for running DQMeasure reads from the ELSE ouput (line 14), rather than the IF output (Line 12).  It seems logical to me that the IF output should be used if the UNITS_FACTOR=1 is true.?

Appreciate any help.  Python scripts included below.

DQMOutputAnalysisHorizontal.py where we specify the units factor

# NAME: ProcDQMOutput.py DQMOutputAnalysis.py
# PURPOSE: Model boresight errors in Lidar data swaths
# Developer: Ajit Sampath, Remote Sensing technologies Project, SGT, Contracter to USGS EROS
# PROJECT: RST, for ASPRS Guidelines on Geometric accuracy of Lidar data
# Libraries: Standard python libraries plus numpy, scipy, matplotlib, liblas
################################################################################
import sys
sys.path.append(r'C:\Program Files (x86)\Python27-14\Lib\site-packages')#(r'C:\Python27\ArcGIS10.2\Lib\site-packages')
sys.path.append(r'C:\Python27\ArcGIS10.3\Lib\site-packages')


import scipy
import os
import numpy,random,time,math
from scipy import stats
from DQMAnalysisFun_Horizontal import *


##############################################################
## NEXT THREE LINES REQUIRE USER INPUTS
UNITS_FACTOR=1 # Change to .3048 for feet and 1 for meter
DataDirectoryPath=r'E:\ALS\P18-037_2016ALS\6aFlightlines_Class6'
OutputDirectorypath=r'E:\ALS\P18-037_2016ALS\7gDQM_py_Class6'


## Next Few lines are DQM calls and Analysis Steps


fCreateBatchFile(DataDirectoryPath, OutputDirectorypath,UNITS_FACTOR)
os.chdir(OutputDirectorypath)
print 'Finished DQManage'
fGenerateAnalysis(OutputDirectorypath, UNITS_FACTOR,250)

DQMAnalysisFun_Horizontal.py.  Although my units factor is 1, the ELSE output is used.  I can adjust to suit the results I want from DQMeasure, however I am thinking there are other analysis consequences of the script miss-understanding the unit factor.

import scipy
import os
import numpy,random,time,math
from scipy import stats




def fCreateBatchFile(DirectoryPath, OutputDirectorypath,UNITS_FACTOR):
    os.chdir(OutputDirectorypath)
    OFile= open('runDQM.bat','w+')
    OFile.write('set path= %Path%;C:\DQMLidar')
    if (UNITS_FACTOR!=1):
        cm="\nDQManage -i "+ '"'+DirectoryPath +'\*.laz'+'"'+' -step 5 -nearest 8 -neighbors_min 8 -neighbors_max 20 -two_ring -eigen_ratio 0.001"'
    else:
        cm="\nDQManage -i "+ '"'+DirectoryPath +'\*.laz'+'"'+' -step 2 -nearest 10 -neighbors_min 8 -neighbors_max 30 -two_ring -eigen_ratio 0.01"'
   
    OFile.write(cm)
    OFile.close()
    os.system('runDQM.bat')
    os.system('del runDQM.bat')


def fGenerateAnalysis(DirectoryPath, UNITS_FACTOR,MIN_PTS):
    textFileList=[]
    for f in os.listdir(DirectoryPath):
        if f.startswith('check') and f.endswith(".txt"):
            textFileList.append(os.path.join(DirectoryPath, f))  


    OFile= open('SummaryStats.txt','w+')
    OFile2= open('PotentialErrorFile.txt','w+')
    OFile3= open('AcceptableFiles.txt','w+')


    OFile2.write('DQM File'+'\t'+'File Name 1'+'\t'+'File Name 2'+'\t')
    OFile2.write('Max Potential Roll Error'+'\t')
    OFile2.write('Number Sample Points'+'\t')
    OFile2.write('Slope'+'\t'+'Intercept'+'\t'+'Dist_Max'+'\t')
    OFile2.write('Mean'+'\t'+'Std'+'\t'+'RMSE'+'\t'+'Discrepancy Angle'+'\t')
    OFile2.write('Num_HS'+'\t'+'Mean HS'+'\t'+'Std HS'+'\t'+'Max D_HS'+'\t')
    OFile2.write('dX'+'\t'+'dY'+'\t'+'dZ'+'\n')


    OFile3.write('DQM File'+'\t'+'File Name 1'+'\t'+'File Name 2'+'\t')
    OFile3.write('Max Potential Roll Error'+'\t')
    OFile3.write('Number Sample Points'+'\t')
    OFile3.write('Slope'+'\t'+'Intercept'+'\t'+'Dist_Max'+'\t')
    OFile3.write('Mean'+'\t'+'Std'+'\t'+'RMSE'+'\t'+'Discrepancy Angle'+'\t')
    OFile3.write('Num_HS'+'\t'+'Mean HS'+'\t'+'Std HS'+'\t'+'Max D_HS'+'\t')
    OFile3.write('dX'+'\t'+'dY'+'\t'+'dZ'+'\n')


    for textFile in textFileList:
        print textFile
        fProcOutPutFile(textFile,OFile,OFile2,OFile3,UNITS_FACTOR,MIN_PTS)
    OFile.close()
    OFile2.close()
    OFile3.close()
    print 'here'




def fGetCenterlineEq(xyz):
    #xyz is an x-y data in most cases for this function
    # all 'planes' should be replaced by line
    xyzm=xyz.mean(axis=0)
    LineEq=numpy.zeros((1,5)) # 1 X 5 vector of planar parameters  
    covMat=numpy.cov((xyz-xyzm).transpose()) # generate covariance matric covMat
    evv=numpy.linalg.eig(covMat)# get eigen vales and eigen vectors of covMat
    indmin=evv[0].argmin()    # gets the index of Minimum eigen value. Gives direction of most variance


    LineEq[0][0]=evv[1][0,indmin] # Median line eq has same parameters as eigen vectors corresponding to min eigen value
    LineEq[0][1]=evv[1][1,indmin] # Median line eq has same parameters as eigen vectors corresponding to min eigen value
   
    a=evv[1][:,indmin][0];b=evv[1][:,indmin][1]; # Median Line eq. aX+bY+c=0 and a^2+b^2=1
   
    # c is determined as below
    c=-(a*xyzm[0]+b*xyzm[1])
    return a,b,c


def fMAD(X,e):
    med=numpy.median(X)
    mad=numpy.median(abs(X-med))
    in1=[x for x in xrange(0,len(X)) if abs(X[x]-med)/mad < e]
    #print in1
    return numpy.array([x for x in xrange(0,len(X)) if abs(X[x]-med)/mad < e])
   
def fProcOutPutFile(fname,OFile,OFile2,OFile3, UNITS_FACTOR,MIN_PTS):    
    fopen=open(fname,'r+');
    riter=1
    Filenames=[]
    for f in fopen:
        if riter<3:
            try:
                #Filenames.append(f.rstrip().split(':')[1])
                Filenames.append(f.rstrip().rpartition('\\')[2])
            except:
                H=0
                   
            #print Filenames
        riter=riter+1
        if riter>3:
            break


    DQM=numpy.zeros((10000,15))
    riter=0
    for f in fopen:
        DQM[riter,:]=numpy.array([float(x) for x in f.rstrip().split()])
        riter=riter+1
        if(riter>=DQM.shape[0]):
            DQM.resize((riter+riter,15))


    fopen.close()
   
    DQM=DQM[DQM[:,0]<>0]
    dX=0
    dY=0
    if(DQM.shape[0]>100):
        dX=max(DQM[:,0])-min(DQM[:,0])
        dY=max(DQM[:,1])-min(DQM[:,1])
    #print dX,dY
   
    if(DQM.shape[0]>100):
        # The next line does:
        # Ratio of eigenvalues should be less than .001
        # DQM[x,7] is the variance of neighborhood size
        # DQM[x,8]/DQM[x,7] is the ratio to enforce that neighborhoods are not very thin or skinny and they are closer to circular neighborhood
        ind_DQM=[x for x in xrange(0,DQM.shape[0]) if DQM[x,9]/(DQM[x,7]+DQM[x,8]+DQM[x,9]) < .001 and DQM[x,7] < 25/(UNITS_FACTOR*UNITS_FACTOR) and DQM[x,8]/DQM[x,7] > .45 ]
        DQM=DQM[ind_DQM,:]
    # Moving origin to Min X,Y,Z for more stable calculation
        DQM[:,0]=DQM[:,0]-DQM[:,0].min()
        DQM[:,1]=DQM[:,1]-DQM[:,1].min()
        DQM[:,2]=DQM[:,2]-DQM[:,2].min()
        if(DQM.shape[0]>MIN_PTS):


            a,b,c=fGetCenterlineEq(DQM[:,0:2]) # Get center line equation
           
            dist_medial=numpy.array([a*x[0]+b*x[1]+c for x in DQM])


            ind_pl=[x for x in xrange(0,len(DQM)) if math.acos(DQM[x,5])*180/math.pi < 5]
           
           
            ind_pl_inliers=fMAD(DQM[ind_pl,6],6)
            ind_pl=[ind_pl[x] for x in ind_pl_inliers]


            ind_Slp=[x for x in xrange(0,len(DQM)) if math.acos(DQM[x,5])*180/math.pi > 15]
           
           
            ind_Slp_inliers=fMAD(DQM[ind_Slp,6],6)
            ind_Slp=[ind_Slp[x] for x in ind_Slp_inliers]
            N=DQM[ind_Slp,3:6]
            Dx=DQM[ind_Slp,6]


            if (len(N)>30):
                dXYZ=numpy.linalg.lstsq(N,Dx,rcond=None)[0] #SPH playing with rcond cuttoff value
            else:
                dXYZ=numpy.zeros((3,1))
           
            if(len(ind_pl)>3):
                slope, intercept, r_value, p_value, std_err = stats.linregress(dist_medial[ind_pl],DQM[ind_pl,6])
                DA=numpy.array([math.atan(DQM[x,13]/abs(dist_medial[x])) for x in ind_pl if dist_medial[x] > 100]) # The Overlap has to be at least 100 units
                if(len(DA)>10):
                    medDA=numpy.median(DA)*180/math.pi
                else:
                    medDA=0


                OFile.write('**********************************'+'\n')
                OFile.write(fname.rstrip().split('\\').pop()+'\t'+Filenames[0]+'\t'+Filenames[1]+'\n')
                err_max=abs(abs(dist_medial[ind_pl]).max()*slope+intercept)
                OFile.write('Max Potential Roll Error::'+str(err_max)+'\n')
                RMSE=numpy.sqrt(numpy.multiply(DQM[ind_pl,6],DQM[ind_pl,6]).mean())


               
               
                   
                if(medDA>0.1 or math.sqrt(DQM[ind_pl,6].mean()**2 + DQM[ind_pl,6].std()**2)>.15/UNITS_FACTOR):
                    print DQM[ind_pl,6].max(),DQM[ind_pl,6].min()
                    OFile2.write(fname.rstrip().split('\\').pop()+'\t'+Filenames[0]+'\t'+Filenames[1]+'\t'+str(err_max)+'\t'+str(len(ind_pl))+'\t')
                    OFile2.write(str(slope)+'\t'+str(intercept)+ '\t'+str(abs(dist_medial[ind_pl]).max())+'\t')
                    OFile2.write(str(DQM[ind_pl,6].mean())+'\t'+str(DQM[ind_pl,6].std())+'\t'+str(RMSE)+'\t'+str(medDA)+'\t')
                    if(len(ind_Slp)>2):
                        OFile2.write(str(len(ind_Slp))+'\t'+str(DQM[ind_Slp,6].mean())+'\t'+str(DQM[ind_Slp,6].std())+'\t'+str(DQM[ind_Slp,6].max())+'\t')
                        OFile2.write(str(dXYZ[0])+'\t'+str(dXYZ[1])+'\t'+str(dXYZ[2])+'\n')
                    else:
                        OFile2.write(str(len(ind_Slp))+'\t'+str('0')+'\t'+str('0')+'\t'+str('0')+'\t')
                        OFile2.write(str(dXYZ[0])+'\t'+str(dXYZ[1])+'\t'+str(dXYZ[2])+'\n')
                       
                           
                else:
                    if(len(ind_Slp)>0):
                        print str(DQM[ind_Slp,13].max())
                    OFile3.write(fname.rstrip().split('\\').pop()+'\t'+Filenames[0]+'\t'+Filenames[1]+'\t'+str(err_max)+'\t'+str(len(ind_pl))+'\t')
                    OFile3.write(str(slope)+'\t'+str(intercept)+ '\t'+str(abs(dist_medial[ind_pl]).max())+'\t')
                    OFile3.write(str(DQM[ind_pl,6].mean())+'\t'+str(DQM[ind_pl,6].std())+'\t'+str(RMSE)+'\t'+str(medDA)+'\t')
                    if(len(ind_Slp)>2):
                        OFile3.write(str(len(ind_Slp))+'\t'+str(DQM[ind_Slp,6].mean())+'\t'+str(DQM[ind_Slp,6].std())+'\t'+str(DQM[ind_Slp,6].max())+'\t')
                        OFile3.write(str(dXYZ[0])+'\t'+str(dXYZ[1])+'\t'+str(dXYZ[2])+'\n')
                    else:
                        OFile3.write(str(len(ind_Slp))+'\t'+str('0')+'\t'+str('0')+'\t'+str('0')+'\t')
                        OFile3.write(str(dXYZ[0])+'\t'+str(dXYZ[1])+'\t'+str(dXYZ[2])+'\n')
    else:
        print 'Not Enuf samples'

[Sam Hackett]

Hi,

An update for records; the author of the scripts has been in contact and clarified that the IF statement in Line 11 is a 'not equal to' statement with the operator != so the ELSE output should be read if the units are equal to 1.  The issue is simply my lack of python knowledge.

Thanks Ajit

Sam