I have a working code for a project, but right now to run through the data it would take years to be completed. Please help me to improve it:

import laspy
import laspy.file
import numpy as np


header = laspy.header.Header()
inFile2 = laspy.file.File("C:\\Users\\Geri\\Desktop\\Sync\\CODETEST\\inverz2.las", mode = "r")
inFile3 = laspy.file.File("C:\\Users\\Geri\\Desktop\\Sync\\CODETEST\\inverz3.las", mode = "r")
inFile = laspy.file.File("C:\\Users\\Geri\\Desktop\\Sync\\CODETEST\\final.las", mode = "rw")


point_records = inFile2.points
point_records = inFile3.points


listx1=np.array([], dtype=float)
listy1=np.array([], dtype=float)
listz1=np.array([], dtype=float)
while t < 1155817:
    while z==0:

        print t
        for points in inFile2.points:
            if ax > xmin and ax < xmax and ay < ymax and ay > ymin:

Till this point the code tries to find the closest neighbor for the input point. The method only works for a data with a few hundred points, but it's too slow for a bigger one.


                                print inFile3.x[t]
                                print inFile2.x[n]
                                newx = [inFile3.x[t]-((inFile3.x[t]-inFile2.x[n])/2)] 
                                newy = [inFile3.y[t]-((inFile3.y[t]-inFile2.y[n])/2)]


                                processz = newy-newx

                                newz = [inFile3.z[t]-((inFile3.z[t]-inFile2.z[n])*ratioz)]


                                listx1=np.append(listx1, (newx))
                                listy1=np.append(listy1, (newy))
                                listz1=np.append(listz1, (newz))

                print n
        if q>0 or p==1:            

inFile.x = listx1
inFile.y = listy1
inFile.z = listz1

Unfortunately I can't install the pyflann or the scipy (for k neighboor) to use them so that's why I tried this workaround. I'm wondering if it's posibble to create a buffer around the selected point and just select by that polygon the new point (like in ArcGIS) (FYI I have an ArcGIS 10.2.2 and a standalone python 2.7)

Description: The code ultimate purpose: (gonna label the code with chapters) get 2 las files (ch1) take one of the point from the first one (ch2), find the closest point from the second file (ch3), create a new point based on the two input points (ch4), calculate new z value adjusted by how far is the new point from the edge of the dataset (ch5), save the new point into a new las file (ch6)

  • 2
    When working with large LAS files, there is going to be more processing time required. That being said, you have nested loops and print statements. If you don't need the print statements (i.e. you already know the process works) you could remove them and that should significantly improve the processing time. You could report at every 10K iterations with something like this: if not t % 10000: print t
    – crmackey
    Oct 2, 2015 at 12:36
  • 3
    My favorite trick for large/very-large processing is adding a periodic check for a "stop.now" file, so that a graceful end to processing can occur it the host needs to be rebooted or the core processing script needs to be updated. Designing the script to be reentrant is a bonus.
    – Vince
    Oct 2, 2015 at 13:25
  • @crmackey Thanks, but while your method makes the process faster, the fact that I run through million points million times makes this code not viable. Right now I'm thinking about to call the Arcpy in and use the Make LAS Dataset Layer (Data Management) ; Select Layer By Location (Data Management) and Describe (arcpy) functions somehow.
    – Gary
    Oct 2, 2015 at 15:13
  • @Gary I think you're on the right track. Maybe a good start would be to create a fishnet for the extent of your LAS dataset and make individual las layers within one cell block of the fishnet at a time (maybe use the extent environment and pass in the extent of each feature of the fishnet as you pass through to build each LAS layer).
    – crmackey
    Oct 2, 2015 at 16:20
  • 1
    In that case I think you should develop a code snippet that focuses on just that critical part and present just that, instead of your whole code, here. See meta.gis.stackexchange.com/questions/4006/…
    – PolyGeo
    Oct 6, 2015 at 23:41


Browse other questions tagged or ask your own question.