3

I have read a LAS file. If I import it in Global Mapper it says the projection, however when I read it with laspy (header.proj_id_1) there is no information associated. I have edited the coordinates and updated the X,Y,Z values to new values (according to this) for the LAS file and now I would like to give a new projection after updating the coordinates. What function do I use?

enter image description here

9
  • 1
    The link scales and offsets, it does not change the projection, it sounds like you're swimming in dangerous waters adjusting the coordinates if you're not 100% sure. If you want to project a las file LASTools is a viable and reasonably cheap option using LAS2LAS which will embed the new CRS in the header. Commented Mar 15, 2021 at 0:23
  • @MichaelStimson I am 100% sure. I use files from the national cadastre agency, that give parameters for the correct coordinates conversion from ETRS89 to Stereo70. The input LAS had UTM zone 35 coordinates, I changed those into ETRS89, then into Stereo70. I have to use the files given by the cadastre agency to make sure the conversion is correct. Other tools, even well known software, do not use the same parameters. The parameters could also change, so I read the parameters every time. I have given the input las the new coordinates, however I need to change the name of the projection. Commented Mar 15, 2021 at 12:46
  • 1
    It seems there are multiple ways to store the CRS according to the spec, so we could use some more information-- what LAS version and point format are you using? Posting a the output from lasinfo for the input and output file might be helpful as well
    – mikewatt
    Commented Mar 18, 2021 at 17:27
  • 1
    Which laspy version are you using?
    – Mapperz
    Commented Mar 19, 2021 at 3:31
  • 1
    I'm not sure how to do it offhand but it is possible-- LAS 1.2 uses the GeoTIFF standard as shown on page 11 here: asprs.org/wp-content/uploads/2010/12/asprs_las_format_v12.pdf But proj_id_1 relates to the "project id", not projection. The VLRs will be the relevant parts of the lasinfo report, I don't see those in your image. It might be a little bit of work to update the georeferencing info, though, since I don't think laspy provides any mechanisms to help with that. Sorry I don't have an easy answer
    – mikewatt
    Commented Mar 23, 2021 at 18:16

2 Answers 2

1
+50

I found a minute to dig into this. We don't quite know what your specific file contains, but here's a quick and dirty way to update the GeoKeyDirectoryTag VLR (page 11) using values as described in the geotiff spec.

This does not currently account for the storage of individual projection parameters in a separate GeoDoubleParamsTag or GeoAsciiParamsTag record. I don't have any example data on hand that contains those params, so if you need help extending this I might need you to post something for me.

This updates the file in place, so be warned.

import laspy

src_path = r"C:\example\test.las"

# see here for possible values:
# http://geotiff.maptools.org/spec/geotiff2.7.html
replace_keys = {
    3072: 26912,  # set ProjectedCSTypeGeoKey to NAD83 / UTM 12N
    3076: 9001,   # set ProjLinearUnitsGeoKey to meters
    4099: 9001    # set VerticalUnitsGeoKey to meters
}


def update_geo_keys(path, lookup):
    with laspy.file.File(path, mode='rw') as f:
        for i, vlr in enumerate(f.header.vlrs):
            # search for the GeoKeyDirectoryTag record
            if vlr.record_id != 34735:
                continue

            modified_body = []
            j = 0

            # loop over geo key entries, stored in groups of 4 values
            while j < len(vlr.parsed_body):
                key_id, tag_loc, count, value = vlr.parsed_body[j:j+4]

                if j > 0 and tag_loc != 0:
                    raise NotImplementedError('modification of projection '
                                              'parameters not implemented')

                # update the value for this key ID, if necessary
                try:
                    value = lookup[key_id]
                except KeyError:
                    pass

                # maintain a new list of values to be written
                modified_body.extend((key_id, tag_loc, count, value))
                j += 4

            # update the vlr with our new values, quit searching
            vlr.parsed_body = modified_body
            break
        else:
            raise ValueError('no GeoKeyDirectoryTag found in file')

        f.header.save_vlrs()


update_geo_keys(src_path, replace_keys)
8
  • Thank you, however I get an error: "NotImplementedError: modification of projection parameters not implemented". Is it somehow possible to delete the projection information, if replacing it is not possible? Commented Mar 25, 2021 at 17:04
  • I found that if I comment the NotImplementedError if statement, it works. Thank you Commented Mar 25, 2021 at 17:20
  • 1
    Glad it's working! The implication of commenting out that exception, however, is that there may be another tag/VLR that has the old projection information still intact which wouldn't match what you've set. That may or may not matter to you. Deleting all projection info would be easy as well, but probably not in-place as it would change the size of the file, requiring it to be rewritten.
    – mikewatt
    Commented Mar 25, 2021 at 17:53
  • I turned your post into a function "def change_proj(src_path):" then the full text, without where you defined the src_path and I get an error, when running it as a function with the path to the las file as argument. error is "[Errno 22] Invalid argument: path/to/file" Commented Mar 29, 2021 at 11:29
  • File "path_to_py_file", line 68, in update_geo_keys f.header.save_vlrs() File "pythonfolderpath\lib\site-packages\laspy\header.py", line 1310, in save_vlrs self.writer.save_vlrs() File "pythonfolderpath\lib\site-packages\laspy\base.py", line 1045, in save_vlrs self.set_vlrs(self.vlrs) File "pythonfolderpath\lib\site-packages\laspy\base.py", line 1068, in set_vlrs self.data_provider.open("w+b") File "pythonfolderpath\lib\site-packages\laspy\base.py", line 144, in open self.fileref = open(self.filename, mode) OSError: [Errno 22] Invalid argument: 'in_path' Commented Mar 29, 2021 at 11:58
0

It sounds like you know the projection, it's just not being read in or possibly stored in the las file correctly? If you want to set the known CRS using laspy, you can like so:

from pyproj import CRS
import laspy

# add known crs
las = laspy.read("original_file.las")
crs_pyproj = CRS.from_string("epsg:4326")
las.header.add_crs(crs_pyproj)
las.write("las_with_crs_added.las")  # save file with crs added to header

# Once set, you can read it back with the built-in laspy method `.parse_crs()`.
laspy.read("las_with_crs_added.las")
epsg = las.header.parse_crs().to_epsg()
print(epsg)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.