coordinate systemellipsoidgeoidinterpolationlidar
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.
How can this interpolation be done?
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.interp2d.html
You have two tasks: (a) find the grid data for the vertices of a grid cell in which a GPS location (x,y) lies and (b) interpolate those data.
An efficient way to accomplish (a) is to store the grid as an array with a single index. Often grids are stored by rows with data progressing left to right within each row. A grid is determined by the coordinates of its origin, (x0,y0) (which will be (-180, -90) for a worldwide grid in decimal degrees), its cellsize c (which is 1/24 in your case), and the number of points per row (n). The values would be recorded in the sequence
(x0,y0), (x0+c,y0), (x0+2c,y0), ..., (x0+(n-1)*c,y0)
(x0,y0+c), (x0+c,y0+c), (x0+2c,y0+c), ..., (x0+(n-1)*c,y0+c)
...
(x0,y0+(m-1)*c), (x0+c, y0+(m-1)*c), ..., (x0+(n-1)*c,y0+(m-1)*c)
This sequence can be addressed with a single index k = 0, 1, ..., n*m-1.
Given an arbitrary point (x,y), you use simple calculations to find the indexes where its four nearest neighbors are stored: they lie between rows i and i+1 and columns j and j+1 where
i = Floor((y-y0)/c), j = Floor((x-x0)/c).
The corresponding offsets are ni + j, n(i+1) + j, ni + (j+1), and n(i+1) + (j+1).
This assumes x0 <= x < x+nc and y0 <= y < y+mc, which is worth checking in a general-purpose program. The coordinates of these four points are readily seen to be equal to
(x0 + j*c, y0 + i*c), (x0 + j*c, y0 + (i+1*c), etc.
Now that you can solve (a) efficiently, it remains to interpolate the values found, which I will denote z00, z01, z10, and z11, situated as suggested in this crude picture (and shown more precisely in the figure at the end):
z01 z11
z
z00 z10
You wish to estimate the value of z. It is convenient to use bilinear interpolation. The formulas amount to the following: first interpolate linearly along the top and bottom edges of the square. Here, only the x coordinate is changing (from x0+j*c to x0+(j+1)*c). Rescale this to t varying from 0 to 1, where
t = (x - (x0+j*c))/c = (x-x0)/c - j.
The interpolated values are
z0 = (1-t)*z00 + t*z10, z1 = (1-t)*z01 + t*z11
Now interpolate vertically between z0 and z1. First rescale y to s varying from 0 to 1, where
s = (y-y0)/c - i
and find the interpolated value
z = (1-s)*z0 + s*z1.
Let (x0,y0) = (-180,-90), c = 1/24 = 0.041667, n = 24*360 = 8640, and (x,y) = (54.4786674627, 17.0470721369) (that is, near latitude 17 degrees north and longitude 54.5 degrees east). Then
i = Floor((y-y0)/c) = Floor(24*(90 + 17.047...)) = 2569
j = Floor((x-x0)/c) = Floor(24*(180 + 54.478....)) = 5627.
The undulation for (i,j), namely z00, is found at offset 8640*i+j = 48,619,849. The undulation for (i+1,j), namely z01, is found at offset 8640*(i+1)+j = 48,628,489. The undulation for (i,j+1), namely z10, is found at offset 48,619,850, and the undulation for (i+1,j+1), namely z11, is found at offset 48,628,490. Suppose these have the values given; that is,
z00 = 31.945, z01 = 31.993, z10 = 31.866, and z11 = 31.911
To interpolate, compute
t = (x-x0)/c - j = 0.488019
s = (y-y0)/c - i = 0.129731
Therefore
z0 = 31.96842 and z1 = 31.88796
Finally
z = 31.958.
Update:
Depending on the ArcGIS version, some LiDAR tools will allow reprojecting through the Environment settings: Output Coordinate System.
For example, the Tile LAS tool:
The output LAS tiles can be reprojected by ... specifying an Output Coordinate System in the environment settings ... .
Original answer:
LiDAR (.las) files and LASdataset files (.lasd) cannot be reprojected within ArcGIS.
A .lasd file is a LAS dataset file developed by ESRI.
The LAS dataset does not import point data contained in the LAS files; it simply stores reference to these LAS files and surface constraints.
According to ArcGIS Help 10.1 one can set a spatial reference to a LAS dataset going to the ArcCatalog, LAS Dataset properties, and then, XY Coordinate System tab.
However, this won't change/reproject the actual .las files. It is only recommended when the .las files do not originally contain the spatial reference and one knows which CRS they belong to.
Strongly recommend reading: LAS dataset considerations - Spatial reference, which includes ESRI's constraints about this topic.
Alternatively, there are other software which will allow reprojection of .las files (example).
Best Answer
PDAL can do this with its filters.reprojection capability, which is based on the vertical datum transformation capabilities of Proj.4.
Make sure that egm08_25.gtx, defined in
GDAL_DATA'svertcs.csvfile asEPSG:3855is available in/usr/local/share/gdalor wherever your system'sGDAL_DATAis defined.PDAL is installable via OSGeo4W or Docker. See https://www.pdal.io/download.html for details.