Manually reversing the rotation should do the trick; there should be a formula for rotating spherical coordinate systems somewhere, but since I can't find it, here's the derivation ( ' marks the rotated coordinate system; normal geographic coordinates use plain symbols):
First convert the data in the second dataset from spherical (lon', lat') to (x',y',z') using:
x' = cos(lon')*cos(lat')
y' = sin(lon')*cos(lat')
z' = sin(lat')
Then use two rotation matrices to rotate the second coordinate system so that it coincides with the first 'normal' one. We'll be rotating the coordinate axes, so we can use the axis rotation matrices. We need to reverse the sign in the ϑ matrix to match the rotation sense used in the ECMWF definition, which seems to be different from the standard positive direction.
Since we're undoing the rotation described in the definition of the coordinate system, we first rotate by ϑ = -(90 + lat0) = -55 degrees around the y' axis (along the rotated Greenwich meridian) and then by φ = -lon0 = +15 degrees around the z axis):
x ( cos(φ), sin(φ), 0) ( cos(ϑ), 0, sin(ϑ)) (x')
y = (-sin(φ), cos(φ), 0).( 0 , 1, 0 ).(y')
z ( 0 , 0 , 1) ( -sin(ϑ), 0, cos(ϑ)) (z')
Expanded, this becomes:
x = cos(ϑ) cos(φ) x' + sin(φ) y' + sin(ϑ) cos(φ) z'
y = -cos(ϑ) sin(φ) x' + cos(φ) y' - sin(ϑ) sin(φ) z'
z = -sin(ϑ) x' + cos(ϑ) z'
Then convert back to 'normal' (lat,lon) using
lat = arcsin(z)
lon = atan2(y, x)
If you don't have atan2, you can implement it yourself by using atan(y/x) and examining the signs of x and y
Make sure that you convert all angles to radians before using the trigonometric functions, or you'll get weird results; convert back to degrees in the end if that's what you prefer...
Example (rotated sphere coordinates ==> standard geographic coordinates):
southern pole of the rotated CS is (lat0, lon0)
(-90°, *) ==> (-35°, -15°)
prime meridian of the rotated CS is the -15° meridian in geographic (rotated 55° towards north)
(0°, 0°) ==> (55°, -15°)
symmetry requires that both equators intersect at 90°/-90° in the new CS, or 75°/-105° in geographic coordinates
(0°, 90°) ==> (0°, 75°)
(0°, -90°) ==> (0°,-105°)
EDIT: Rewritten the answer thanks to very constructive comment by whuber: the matrices and the expansion are now in sync, using proper signs for the rotation parameters; added reference to the definition of the matrices; removed atan(y/x) from the answer; added examples of conversion.
EDIT 2: It is possible to derive expressions for the same result without explicit tranformation into cartesian space. The x, y, z in the result can be substituted with their corresponding expressions, and the same can be repeated for x', y' and z'. After applying some trigonometric identities, the following single-step expressions emerge:
lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ
It depends on what version of the LAS specification you are using. If it is 1.3 or less, then the specs define georeferencing information using pre-defined (see specs) variable length records (VLRs) using the same format as the GeoTIFF:
Georeferencing for the LAS format will use the same robust mechanism
that was developed for the GeoTIFF standard.
This format, though somewhat challenging to grok at times, is remarkably flexible. It relies on three defined tags called the GeoKeyDirectoryTag tag, which is like a table of contents for georef data, the GeoDoubleParamsTag tag, which is like a store of all double-precision values referred to in the GeoKeyDirectoryTag, and the GeoAsciiParamsTag tag, which similarly is used to store all ASCII (text) values. This site provides a good explanation and an example.
As of LAS v. 1.4 however, this method of storing georeferencing information was changed to favour the well-known text (WKT) format, also stored in defined VLRs, although the GeoTIFF format is still used for legacy:
The Coordinate Reference System (CRS) information for the point data
is required for all data. The CRS information will be placed in
Variable Length Records or Extended Variable Length Records (note that
if the writer wishes to maintain legacy compatibility, then GeoTIFF in
VLRs must be used). The CRS is represented by either GeoTIFF or Well
Know Text as indicated by the WKT Global Encoding bit. Point Record
Formats 0-5 can use either GeoTIFF or WKT (but not both
simultaneously). Point Record Formats 6-10 must use WKT.
I see no reason given these flexible formats why you couldn't store point info in geographic coordinates (lat/long) but this would be fairly unusual for LAS data in that I've never seen it done previously. I imagine the reason is that LiDAR datasets tend to be of rather large scale (small spatial extent) and projected coordinate systems are therefore preferred. It makes calculating the distances between points, which is important for some algorithms (e.g. point classification or filtering), much easier.
Best Answer
Going from a lon-lat coordinates (spheroid) to cartesian coordinates is not as easy as just using a mathematical formula, the formulas for the reprojection are quite big and the solution is not explicit sometime.
Depending on the accuracy required there are standarized procedures to undertake the task. A common one is project a lon-lat data into a Transvers Mercator projecton. The central meridian and the scale factor can be choosen at will in order to match the accuracy required. Here some reading about the subject by Richard Stanaway Link.
Once that is said, your problem could be seen as two problems:
For the first task I recommend you a sound GIS tool like the user friendly QGIS or the proj.4 if you want to code.
Many questions are still to come: