Converting between coordinate systems for LAS file using LASPY

Question

I am trying to convert a .las file from WGS84 to UTM17. I have some drone imagery that was used to generate a 3d point cloud/ .las file, but it is in latitude/ longitude instead of easting/ northing.

I tried using laspy, a great python API for reading las files. I was able to use this to convert the XYZ info from the .las to a csv file. I then used this csv file and a package in R to convert the coordinates. Now I am trying to use laspy to write a new .las file with the correct coordinates.

Here's my code so far:

import numpy   
import laspy 
#pandas used here to convert csv to       array
import pandas
#Import a .las file and name it something
van_taken = laspy.file.File("C:\Users\Mainframe2\PycharmProjects\LAS   Conversion\Coordinate Conversion\Taken.Lands.Points.las", mode="rw")
converted_csv = pandas.read_csv("C:\Users\Mainframe2\PycharmProjects\LAS    Conversion\Coordinate Conversion\Taken.Lands.UTM.csv")
converted_csv = converted_csv.as_matrix()
#get individual column from array using Pandas
converted_x = converted_csv[:,1]
converted_y = converted_csv[:,2]
converted_z = converted_csv[:,3]
#Convert original points to new points
van_taken.x = converted_x
van_taken.y = converted_y
van_taken.z = converted_z

This is where the problem arises.

van_taken.x
array([-83.0830127, -83.0830128, -83.0830127, ..., -83.0830927,
   -83.0830926, -83.0830929])
converted_x
array([ 269873.21570411,  269872.62699301,  269873.32416573, ...,
    269073.15716998,  269074.64298861,  269071.18083953])
van_taken.x = converted_x

But then!!!!

van_taken.x
array([-297.8583648, -297.8583648, -297.8583648, ..., -297.8583648,
   -297.8583648, -297.8583648])

For some reason, it appears that when I try to set them equal, van_taken.x refuses to be equal to what it says.

I've looked through the documentation and can't seem to find anything. Is there some scale or offset that is causing this to happen? Once I figure this out I will be able to save the .las file and convert to the correct coordinate system.

>>> van_taken.header.scale

[1e-07, 1e-07, 0.01]
>>> van_taken.header.offset

[-83.11, 10.0, -1000.0]

I changed the scale to

van_taken.header.scale = [1,1,1]

>>> van_taken.x
array([ 269872.89,  269872.89,  269872.89, ...,  269072.89,  269074.89,
    269070.89])
>>> converted_x
array([ 269873.21570411,  269872.62699301,  269873.32416573, ...,

Changing the scale seemed to make the number more realistic, but they are still different now. I'm going to experiment with a couple different scales.

Update

Grant Brown's answer works like a charm. I copied and pasted it. For some reason when I ran the whole code at once in Pycharm, it gave me some syntax errors. But I just ran it line by line and it didn't have a problem.

The results:

>>> converted_x
array([ 269873.2157 ,  269872.62699,  269873.32417, ...,  269073.15717,
    269074.64299,  269071.18084])
>>> van_taken.x
array([ 269873.2157    ,  269872.62699011,  269873.32417004, ...,
    269073.15716974,  269074.64298997,  269071.18084039])

Grant's answer is accurate at least to the number of decimal places we have in our file. (I reduced the number of decimals I originally had from 10 places to 5, since 5 is sub milimeter precision.)

Finally, I just want to say remember to close the file. I was a complete noob here, but after closing the file I noticed the file size was exactly the same size. I therefore incorrectly concluded that the file was not overwritten. But, as Grant explains.

"Unless you changed the number of points, the file is guaranteed to remain the same size. It's a binary file, and everything in it takes up a very precise amount of space. No matter what you do, LAS files store coordinates as integers. This is not something you can change or convert, ever. What you can do is change the scale and offset value, which allows you to construct coordinates in whatever system you want; they're the map between the data storage format and the coordinates you're actually interested in."

I got this one sorted out.

Answer 1

You should check what the scale and offset are for your file. This can be done as follows:

van_taken.header.scale
van_taken.header.offset

This almost looks like an overflow error to me. The lower case x, y, and z properties need to re-scale and re-offset the coordinates to store it as an integer (which is how LAS files store them). To be honest, setting the values of the coordinates as scaled values is a bit of an anti-patern, because you lose control of how precision is lost.

Here's the code that handles your set operation, located in base.py:

    def set_x(self, X, scale = False):
        '''Wrapper for set_dimension("X", new_dimension)'''
        if not scale:
            self.set_dimension("X", X)
            return
        self.set_dimension("X", np.round((X -     self.header.offset[0])/self.header.scale[0]))
        return

If you can figure out what the unscaled integer representation of your converted data should be, that's probably a better way to store it (e.g., using the capital letter X, Y, and Z properties of the file).

If you're fine with the above approach to converting between the integer representation and the floating point representation, then I'd consider adjusting your scale to ensure that you don't end up with integers greater than four bytes in size.

If this isn't explicable via integer overflow due to scaling issues, then we definitely need to figure out what's going on. If it's an overflow issue, I'd be open to trying to guard against this case so long as it doesn't have too terrible a performance penalty.

Edit:

It looks like overflow is definitely the issue.

When you're assigning a scaled coordinate value into a LAS file, laspy needs to find some way of representing it as a four byte integer. Currently, it faithfully believes the information in the header. That is, it will subtract the appropriate offset (for the X, Y,or Z dimension) and divide by the appropriate scale (for the X, Y, or Z dimension). The result is then rounded to produce an integer.

Your file has an X scale of 1e-7, and an X offset of -83.11. Thus, to convert any new scaled value of x to its integer representation (which is what happens when you assign into the lower case 'x' property of your file), you need to add 83.11 and divide by 1e-7. For your first value, 269873.21570411, this results in a value of 2.699563e+12. The largest number you can store in four bytes is 2.14e9 for signed integers and 4.29e9 for unsigned.

Currently, laspy doesn't check for this mistake, resulting in an integer overflow. As I mentioned above, it's probably best to assign the integer values (to the capital X, Y, and Z properties) yourself to avoid any ambiguity.

As a quick fix, however, you can simply change the offset. The following ought to work:

van_taken.header.scale = [0.01,0.01,0.01]
van_taken.header.offset = [0,0,0]

You can increase the precision of your conversion by using a large offset and large scale. For example, if all of your scaled X coordinates are greater than 200000, you could use an X offset of 200000. Then, when a small scale like 1e-7 is used, the numbers it will be inflating will be smaller. That's something to play around with, keeping in mind the four byte limit.

In a lot of problems, and in a lot of computing environments, it's easy to gloss over the fact that floating point arithmetic is fundamentally not like real number arithmetic. Unfortunately, working with LAS files is not one of those cases.

Edit 2:

So can I change the scale value? Will this affect the data in any other ways?

Yes, the reason you can change the scale in this case is that you're supplying the scaled value. If you tell the LAS file that the scale is some particular value, that's the value it will use when re-scaling data. You wouldn't want to mess with the scale if you were reading existing LAS data.

Edit 3

Last question. So it seems that you were right about the scales. But I tried the .01 and it wasn't giving me the number. I then tried 1.000001 and it seemed closer, but they still aren't the same. Any tips for selecting the correct value, other than trial and error?

I don't think there's really a 'correct' value per-se, it's a trade off. No matter what you do, you're trying to store a (probably 8 byte) floating point vector as a vector of 4-byte integers multiplied by an 8 byte (double precision) floating point scale term and added to an 8 byte floating point offset term. That's not a lossless conversion in general, but it's what you have to do to store data in a standards compliant LAS file.

I would consider trying something like this (not tested):

converted_x = converted_csv[:,1]
converted_y = converted_csv[:,2]
converted_z = converted_csv[:,3]
xmin = np.floor(np.min(converted_x))
ymin = np.floor(np.min(converted_y))
zmin = np.floor(np.min(converted_z))

xmax = np.ceil(np.max(converted_x))
ymax = np.ceil(np.max(converted_y))
zmax = np.ceil(np.max(converted_z))

xrg = xmax - xmin
yrg = ymax - ymin
zrg = zmax - zmin

van_taken.header.offset = [xmin,ymin,zmin]
safety_factor = 2
maxval = 2e9
van_taken.header.scale = [safety_factor*xrg/maxval, safety_factor*yrg/maxval, safety_factor*zrg/maxval]    
van_taken.x = converted_x
van_taken.y = converted_y
van_taken.z = converted_z