8

Is there any tool to convert all Z attributes (ellipsoid h) value into orthometric height (H)? The idea is, having a point cloud (.las) and geoid undulation EGM08 (.asc) get a point cloud (.las) with Z orthometric values.

To get orthometric heights is necessary calculate N (geoid height) and is required using a bi-linear interpolation.

enter image description here

How can this interpolation be done?

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.interp2d.html

Improve this question
10
  • Terrascan can do it, though I'm not sure how to push all the buttons to do that. Where did you get the LAS from? Can you ask your supplier to give the LAS files to you in orthometric height? Isn't it EGM08 (Earth Gravitational Model) not EMG08? It may seem like being pedantic but it could make a big difference in searching for the right term.
    – Michael Stimson
    Commented Oct 8, 2017 at 21:10
  • 1
    There's lots of ways to attack this problem. Assuming you want to stay in Python, if you have your own grid files, you could read up the points with LASpy, read the grid file using GDAL, and write out a new file with LASpy representing your orthometric heights. If you want to simply run a command, PDAL can convert ellipsoid (EPSG:4326) LAS data to orthometric (EPSG:5773) using the vertical datum transformation capabilities of Proj.4. The question needs more detail on what you have and what you are trying to achieve.
    – Howard Butler
    Commented Oct 9, 2017 at 14:22
  • Have you looked into VDatum from NOAA? vdatum.noaa.gov
    – Aaron ♦
    Commented Oct 9, 2017 at 15:15
  • I have a .las from LiDAR with more than 4M points. I check the height of these points comparing with a surface model, and I find out that Z values are elipsoid height. I need to transform Z values into a orthometric. @MichaelStimson you are in the path, I made a python script that reads a .asc (EGM08) and construct a function ((scipy.interpolate.interp2d) to get a geoid height giving latitude and longitud. I checked the results with a few points and it's ok. To avoid make a module for this task, I asked for a tool that can make this. But I think I am very close to the solution with LASpy.
    – BenAC
    Commented Oct 9, 2017 at 19:27
  • Fair enough. Although Terrascan can do it I wrote my own C++ console program to do what you're doing in python; 4 million points (or was that pulses?) is a very small number of records when it comes to LAS, so you should be able to process that in python but if this is something you're doing regularly or for larger areas you should consider some LiDAR processing software like Terrascan, mostly for fixing when you find anomalies in the LiDAR that affects your output.
    – Michael Stimson
    Commented Oct 9, 2017 at 21:21

1 Answer 1

Reset to default
6

PDAL can do this with its filters.reprojection capability, which is based on the vertical datum transformation capabilities of Proj.4.

pdal translate input.las output.las reprojection 
    --filters.reprojection.in_srs="EPSG:4326+4326" 
    --filters.reprojection.out_srs="EPSG:4326+3855"

Make sure that egm08_25.gtx, defined in GDAL_DATA's vertcs.csv file as EPSG:3855 is available in /usr/local/share/gdal or wherever your system's GDAL_DATA is defined.

PDAL is installable via OSGeo4W or Docker. See https://www.pdal.io/download.html for details.

Improve this answer
3
  • 1
    Thanks first of all. What is *.gtx? I´ve never seen this before. I use a .asc .asc, may be is usefull for this problem because y use it and the resoults are perfect. In this resource I can´t see this *.gtx format. In this pdal function, the egm grid is to interpolate is fixed in vertex? Where can you define the interpolation method?
    – BenAC
    Commented Oct 10, 2017 at 9:34
  • GTX is the raster format that Proj.4 uses for grid shift files.
    – Howard Butler
    Commented Oct 10, 2017 at 14:58
  • FWIW, I have packaged PDAL 1.6.0 including vertical datums here: copr.fedorainfracloud.org/coprs/neteler/pdal (I hope I got it right :-)
    – markusN
    Commented Oct 29, 2017 at 10:41

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.