FOSS4G 2018: wait, what?

In late August I went to FOSS4G2018, held in conjunction with the HOTOSM summit in Dar Es Salaam, Tanzania. I’ve been privileged to attend three of these events in a row now, and this one is incredibly difficult to wrap up neatly and put a bow around.

Why is this so?

Phew. The first question is where to even start? To begin with, it was held in the city of Dar Es Salaam, on the east coast of Africa. Flying in on my comfortable 777-300 I literally had no idea what to expect – it was my first visit to the cradle of humanity.

And wow. Straight off the plane, Africa launched itself with an exuberance I was completely unprepared for.

…which set the scene for the week ahead. I’d arrived early, because I was presenting a workshop about PDAL on the first day of the week. This was a story in itself, with my co-presenter blocked from travel literally at the last minute!

I’m digressing. Back to linear time storytelling. I felt completely overwhelmed by Africa to start with, and then – that state continued. So let’s wrap some boundaries around it using days.

Sunday.
I had Sunday to get my bearings both in Dar Es Salaam and the PDAL workshop. It’s the first time I’ve given it, so the usual nerves were well and truly entrenched. I set off toward the conference centre from my hotel; and was greeted by a jovial rastafarian who named himself Picasso. I wandered around with him for more than an hour, learning about Dar, about Tanzania, about how people live. He was joined by a friend of his, and at some point, I knew it was going to turn into a transaction.

Sure enough, some paintings appeared and it was time to do business. I bought one – it seemed really good value, we’d had a great chat about all kinds of things. My Australian conservativeness says ‘ah you were hustled good!’; but my human inner core says ‘well, that was meant to be – you learned a lot in that time, and connected with other humans from someplace completely foreign, and you helped each other out in ways which are both immediate and not so immediate’.

Hakuna matata, says the universe.

So I returned to my wandering, found the conference centre then went back to prepare my workshop.

…which went well, making sure everything mostly worked right up until the post-Geochicas icebreaker event, which happened late-ish at night in a rooftop bar – the to-become-infamous High Spirits lounge. First, I urge you to read this article about the importance of the Geochicas women’s meetup. Next, to the men who turned up – please note how it felt to take that step and respond to the invitation. Often, by default, tech and other events are ‘male focussed’ – so if we felt somewhat overwhelmed, remember it. It’s a reminder of how minority groups can feel all the time.

So I stayed up too late (of course) meeting the Geochicas and friends.

The workshop

Monday bought the PDAL workshop. With 20 people in the room I was terrified I’d do the material no justice at all – especially given a last minute halving of the facilitating team. But we made stuff, people learned things, I have a lot of feedback to work with, and can help improve future editions of the workshop material and my own workshop presentation skills.

One of the awesome things about running workshops is that the room is filled with people who have different problems to the ones I work on all the time. It’s an amazing exercise to both deliver information which is hopefully useful; and also come across challenges ‘in flight’ to try and work with. I was incredibly flattered to see people applying their new knowledge more or less straight away!

If you’ve ever considered a workshop and thought ‘well, I’m not ready yet..’ – then do it! Your use case is your own, you use tools in a way few if any other people do, and your experience is extremely valuable to the community. You’ll also learn new tricks in the process.

The rest of my day was spent in an OpenDroneMap workshop; and later that evening, checking out the Badminton society for an informal social dinner.

This was an experience in itself – it actually was a badminton club! But we could sit out in a courtyard and order local beers with incredible indian food. The night passed, conversations were had, walking buses took people home. #community

Tuesday. A day off.. or not.

Tuesday – I had planned a day off but went to see a HOTOSM discussion about issues faced by women mapping in Africa. It was a fantastic discussion, great to sit in and listen to. As I wrote notes, I realised many of the issues faced by women in Africa also occur in Australia. However, for African women the threat of physical violence is far more real; with fewer support mechanisms – especially further from cities.

I learned about women mapping for basic infrastructure – for example where safe places are; or where people can go for help. These essential services, things which create a civil society; are being built at a grassroots level. Amazing stuff! To repeat a link, I urge you to read this article, which recaps a lot of the discussion – what can we do to better support female mappers. I think these can apply in any country!

Afterward, a small crew of tall white people descended on the Kariakoo markets. It was intense! A bustling street market precinct full of life and people and everything. I bought a coconut (nearly a ritual by now.. the daily coconut) before we escaped, completely overwhelmed.

This was followed by an icebreaker function at … the Badminton society. So more curry and beer and amazing conversations. Late. Into. The. Night… and walking buses for whoever needed them.

Wednesday – and the conference proper!
By this point I was much relieved, my job done and now time to settle in and listen, learn, and talk to people. I was also feeling OK about walking around Dar Es Salaam – I’d opted to travel by foot every day, roughly 20 minutes walk through the city centre. On foot, the city is quite incredible – it opens gently like a flower in the morning, hums madly all day, and gently settles into a long evening.. with people out on the streets talking, mixing, eating until well into the night. An interesting hustle, but in a rhythm; and a purposeful pace.

So the conference. A most amazing opening address from the Tanzanian minister responsible for the Environment, January Makemba. Pushing heavily into the ‘why’ of mapping; and who we are as an open source community. We do it to make the world better, and we are socially conscious. I agree! I hope you do as well!

For the rest of the day, Eric Lemoine educated me about managing point clouds in Postgres-pointcloud and the LOPoCS/itowns stack; then Steven Feldman provided an insightful reflection on open communities and how we can (always) strive for better. Next, Peiro Toffanin and Dakota Benjamin showed their amazing work on OpenDroneMap and sharing results to OpenAerialMap. I finished the day with another OpenDroneMap presentation – Tomas Holderness deploying ODM on an AWS stack using a simple Python SDK.

At a break, I went out and bought a coconut to drink.

My mind thoroughly overwhelmed, the conference dinner was up next. And then, there was more mind blowingness. We were herded into buses, and whisked off to the dinner venue with a police escort! Once there, we talked, saw some incredible dance and drumming performances, Astrid Emde won a Sol Katz award and then we danced.

Gala dinner views

…finally, we were bussed back to Dar Es Salaam, and kicked back with beer and conversations. Late. Into. The. Night.

Thursday – peak overwhelm. Peak what? You can’t peak yet…
Thursday started with the OSGeo AGM – introducing the first Oceania-based FOSS4G event since 2009!

…then a ‘pub quiz’ – which actually let us explore some of our biases about the world and how it is represented.

I chaired a session with some talks I would never have prioritised – which is always a fantastic exercise in serendipity. The talks were excellent – reconstruction of historical buildings in their geospatial context; open source routing engines; gRPC and protocol buffers for geospatial service building; and machine learning for image recognition.

Then back to the plenary room for a keynote delivering an overview of the Zanzibar Mapping Initiative – which you should definitely check out.

And of course, a daily coconut.

B2B rooftop bar moar dancing and networking.. Late. Into. The. Night.

Friday… getting there…
I saw Tom Van Tillburg create city scale models from PostGIS and PDAL; and peer-to-peer pointcloud sharing systems from Thibault Mutabazi. I learned about an application for ‘cycling mentors’;  the application of open source geospatial software in the World Food Program’s data infrastructure; and how the Copernicus program is integrating open source software in it’s infrastructure. Finishing the day, and the event, were two thoughtful and powerful keynote talks on what diversity, inclusion and community mean from Angela Oduor Lungati of Ushahidi, and Maria Arias de Reyna –  OSGeo president.

Sometime that day, I found a coconut seller and bought a coconut.

And once more, a bit of dancing then discussions Late. Into. The. Night.

*boom*. By now, my mind was properly blown. Time for a holiday?

Zanzibar and recovery/unpacking/days off.

Dhows in the distance.. from the Stonetown patio

…to Zanzibar. the perfect place for a well earned rest, right?

Yep. but nope.

My plan for three days of no internet, yoga, and deep relaxation on Zanzibar were broken. They started as they meant to continue with a quick trip to the OSgeo code sprint and a mad dash to the ferry by tuktuk, ably assisted by Nils Hempelmann. I’d also invited part of the ODM team to stay at the airBnB while they needed to. So we talked shop, ate, drank coffee, drank spiced tea, ate more curries, went to the night markets, concocted more ideas, thought about dancing, slept for more hours than the past three days combined, went to a beach, swam in the ocean, drank beer, talked more shop, became locals at a seafront café with amazing masala tea, went on a spice tour, talked shop some more, tried to organise helping out some drone workshop folks with accomodation, danced, drank more beer, and finally, finally, finally squelched off to the Zanzibar airport and the way home.

I don’t really apologise for the density of that text block – that’s exactly what it felt like.

In fact describing the whole week would have fit way better with a massive blob of unpunctuated text! Like a firehose of consciousness, being directly injected…

By now – my train was thoroughly wrecked and even ten hours in Dubai airport working on the FOSS4G SotM Oceania conference seemed sane.

Do not do this – do not spend ten hours in an airport working. It is a batshit crazy idea! 

Be clever and book your flights more carefully.

In conclusion
I have probably left you with an impression that FOSS4G conferences are pretty insane, and so are the people that go to them (mea culpa).

What I hoped to transmit is that this event is actually an incredible week of learning and communication. For a geospatial practitioner, or researcher, or coder, or even anyone just interested it is an amazing cornucopia of knowledge and experience.

There is intense technical discussion, there is intense reflection on the community and it’s makeup and actions.

It is also intensely social. I am the wrong person to comment on inclusivity and social happenings at these events, because I see through an extremely privileged lens. I can say that I do highly appreciate the difference between, say, socialising at FOSS4G and socialising at, say, any other tech focussed event. There is a tangible difference between your average industry spatial thing and FOSS4G, which January Makamba  summed up well in his introduction:

look at this room – it’s beautiful!

…to a room full of vibrant colours, genders, cultures and ideas. The Dar Es Salaam committee worked on the mission ‘leave no one behind’ – and they did an amazing job!

It is not a perfect community – however, it is prepared to reflect on itself; and act on those reflections. There is support and respect, people who go out of their way to make walking buses so others feel safe to get to their accomodation, people who volunteer their time and energy to create a better world. FOSS4G is a week – one week in the year – where you get to catch up on half conversations that go on all year, to find new connections and strengthen old ones. To learn, to give, to receive, to grow. As a professional and as a human. It’s undeniably an incredible privilege.

FOSS4G in 2018 stood out because of how far it went to encapsulate all of these things.

…and it’s not over yet. You can see slides from presentations; and session videos where available – opening and closing sessions at the time I am writing this. I hope these collections grow!

Asante Sana, Dar Es Salaam; and all the people who make this community; and all the wonderful people I spent hours discussing ideas far beyond geospatial software with.

Finally, a massive thanks to Synthesis Technologies for letting me spend a large chunk of a small travel budget to get there!

Here’s hoping I can make it to Romania for the 2019 edition…

16 bit to 8 bit RGB colours with PDAL and Python

Sometimes LAS files have RGB values stored as 16 bit colours. This currently upsets the Potree viewer, and maybe some other things. The following recipe using Python and PDAL ‘eightbitifies’ colours. It also compresses incoming LAS files to LAZ.

PDAL has a neat filters.python option, which passes every incoming point into a python function. Flexible huh?

First, a Python function to scale colours:

import numpy as np

def eightbitify(colour):
    notzero = np.where(colour > 0)
    colour[notzero] = (colour[notzero]/255) - 1
    return colour

def scale_colour(ins,outs):
    outs['Red'] = eightbitify(ins['Red'])
    outs['Green'] = eightbitify(ins['Green'])
    outs['Blue'] = eightbitify(ins['Blue'])
    return True

ins is a numpy array of incoming points from PDALs reader. PDAL dimensions define what’s in there – so here I’ve asked filters.python to read Red, Green, and Blue into numpy arrays to work on. The entire set of data will be loaded up, making python array operations useful and fast (if you have the memory). PDAL’s —stream operator isn’t enabled for filters.python –  if you are memory constrained, think about filters.splitter to make manageable chunks first.  This actually makes sense – PDAL has no way of knowing what you’re going to write into your python function!

Next, construct a JSON pipeline to run the conversion:

 

{
    "pipeline": [
        {
            "type" : "readers.las",
            "filename" : "file.las"
        },
        {
            "type" : "filters.python",
            "script": "/opt/data/scalecolour.py",
            "function": "scale_colour",
            "module": "scalecolour"
        },
        {
            "type" : "writers.las",
            "filename" : "outfile.laz"
        }
    ]
}

 

…and finally invoke PDAL to make it go. Here it’s wrapped in a shell script to process a bunch of data. Highlighted lines show the actual PDAL invocation. I’ll replace the shell script to loop over files with a python function once it’s written.

for f in lasfiles/*.las;
    do
    #echo $f
    fileout=$(basename $f ".las")
    #echo $fileout
    docker run -it -v /set/your/working/path:/opt/data \
                    pdal/pdal pdal \
                    pipeline /opt/data/scale_colour.json \
                    --readers.las.filename=/opt/data/${f} \
                    --writers.las.filename=/opt/data/lasfiles_8bit/${fileout}.laz

    done

 

I ran this using the pdal/pdal:latest docker image:

docker pull pdal/pdal

 

EGU 2017

I’m going to Vienna next week for EGU17 – which is awesome, and I’m extraordinarily excited under my ‘I’m a badass science ninja’ veneer! I’ll be listening a lot, and talking a little bit about work I’m doing at the National Computational Infrastructure on some ideas for data services around massive, dense point data – standing in front of the poster pictured above.

Hope to see you there!

ps – I promise, I’ll write about sea ice real soon. I need to – for Pint of Science (Australia)  next month…

The LiDAR uncertainty budget part III: data requirements

Part I of this series described the LiDAR georeferencing process (for one, 2D, single return scanner). Part II showed how point uncertainties are computed and why they are useful.

What hasn’t been covered so far is what data you need to do this stuff! Here is your shopping list:

Scanner data

  • Raw LiDAR ranges
  • Scanner mirror angles
  • Waveforms with times if you are going to decide range-to-returns for yourself
  • Intensity
  • anything else (RGB, extra bands, etc)

Essentially you need the full data package from the LiDAR including any corrections for mirror wobble or other sensor oddness – and a way to decode it. You also need engineering  drawings which show the LiDAR instrument’s reference point.

Aircraft trajectory data

  • GPS/GNSS based aircraft trajectory (positions in space)
  • Attitude data  (relationship between the aircraft’s orientation and some reference frame) normally collected by a ‘strapdown navigator’ – which is an inertial motion sensor (often called IMU), and often integrated with GNSS.

Hopefully you don’t have to do the work combining these yourself, but you need a high-temporal-resolution data stream of position (XYZ) and attitude (heading pitch roll, or omega phi kappa, or a quaternion describing rotation of the airframe relative to the earth). This needs to be as accurate as possible (post-processed dual-frequency GPS with preferably tightly-coupled processing of IMU and GPS observations).

Navigation instrument and aircraft data

  • Engineering drawings which show the navigation instrument’s reference point
  • The lever arm between the navigation instrument reference point and the LIDAR reference point, in IMU-frame-coordinates
  • The rotation matrix between the IMU and the LIDAR coordinate systems
  • If necessary, the rotation matrix describing the difference between the aircraft’s body frame and the IMU’s internal coordinate system (quite often, this gets simplified by assuming that the IMU frame and the aircraft frame are equivalent – but many operators still account for this – which is awesome!)
  • Any boresight misalignment information you can get (tiny angles representing the difference between how we think instruments were mounted and how they actually were mounted)

Of course, you could also just push your contractors to deliver point-by-point uncertainty or some per-point quality factor as part of the end product.

This series might have one more post, on how to deal with these data – or more accurately a tutorial on how not to. And then I run out of expertise – brain = dumped, and over to people with newer, more up to date experience.

The LiDAR uncertainty budget II: computing uncertainties

In part 1, we looked at one way that a LIDAR point is created. Just to recap, we have 14 parameters (for the 2D scanner used in this example) each with their own uncertainty. Now, we work out how to determine the geolocation uncertainty of our points.

First, let’s talk about what those uncertainties are. An obvious target is GPS positioning uncertainty, and this is the source that comes to mind first. Using a dual-frequency GPS gives the advantage of sub-decimetre positioning, after some post-processing is done. Merging those positions with observations from the IMU constrains those positions even further.

For a quick mental picture of how this works, consider that the navigation unit is collecting a GPS fix every half a second. The accelerometers which measure pitch, roll and heading take a sample every 1/250th of a second. They also keep track of how far they think they’ve moved since the last GPS fix – plus the rate of motion is tracked, putting a limit on how far it is possible to move between each GPS fix. The navigation unit compares the two data sources. If the GPS fix is wildly unexpected, a correction is added to GPS positions and we carry on.

So we get pretty accurate positions (~5cm).

But is that the uncertainty of our point positions? Nope. There’s more. The laser instrument comes with specifications about how accurate laser ranging is, and how accurate it’s angular encoder is.

Even more, the navigation device has specifications about how accurate it’s accelerometers are, and all of these uncertainties contribute! How?

Variance-covariance propagation to the rescue

Glennie [1] and Schaer [2] used Variance-covariance propagation to estimate the uncertainty in geolocation of LiDAR points. This sounds wildly complex, but at it’s root is a simple idea:

uncert(thing1 + thing2) = \sqrt{uncert(thing1)^2 + uncert(thing2)^2}

See that? we figured out the uncertainty of a thing made from adding thing1 and thing2, by knowing what the uncertainties of things1 and thing2 are.

Going from simple addition to a few matrix multiplications was a quantum leap, handled neatly by symbolic math toolboxes – but only after I nearly quit my PhD trying to write out the uncertainty equation by hand.

Here is the uncertainty equation for LiDAR points, whereU is the uncertainty:

U = FC_uF^t

What!! Is that all? I’m joking, right?

Nearly. F is a 14-element vector containing partial derivatives of each the parameters that go into the LiDAR georeferencing equation (the so-called Jacobians). C is a 14 x 14 matrix on which the uncertainties associated with each element occupy the diagonal.

Writing this set of equations out by hand was understandably a nightmare, but more clever folks than myself have achieved it – and while I certainly learned a lot about linear algebra in this process, I took advice from an old astronomer and used a symbolic maths toolbox to derive the Jacobians.

…which meant that the work actually got done, and I could move on with research!

Now my brain is fully exploded – what do uncertainties look like?

Glennie [1] and Schaer [2] both report that at some point, angular motion uncertainty overtakes GPS position uncertainty as the primary source of doubt about where a point is. Fortunately I found the same thing. Given the inherent noise of the system I was using, this occurred pretty quickly. In the Kalman filter which integrates GPS and IMU observations, a jittery IMU is assigned a higher uncertainty at each epoch. This makes sense, but also means that angular uncertainties need to be minimised (for example, by flying instruments in a very quiet aircraft or an electric-powered UAS)

I made the following map years ago to check that I was getting georeferencing right, and also getting uncertainty estimates working properly. 

It could be prettier, but you see how the components all behave – across-track and ‘up’ uncertainties are dominated by the angular component not far off nadir. Along-track uncertainties are more consistent across-track, because the angular measurement components (aircraft pitch and a bit of yaw) are less variable.

The sample below shows LiDAR point elevation uncertainties (relative to an ITRF08 ellipsoid) during level flight over sea ice. At instrument nadir, height uncertainty is more or less equivalent to instrument positioning uncertainty. Increasing uncertainties toward swath edges are a function of angular measurement uncertainty.

But normally LiDAR surveys come with an averaged accuracy level. Why bother?

i. why I bothered:

In a commercial survey the accuracy figure is determined mostly by comparing LiDAR points with ground control data – and the more ground control there is, the better (you have more comparison points and can make a stronger estimate of the survey accuracy).

Over sea ice this is impossible. I was also using these points as input for an empirical model which attempts to estimate sea-ice thickness. As such, I needed to also propagate uncertainties from my LiDAR points through the model to sea-ice thickness estimates.

In other words, I didn’t want to guess what the uncertainties in my thickness estimates were. As far as plausible, I need to know what the input uncertainty is for each thickness estimate is – so every single point. It’s nonsensical to suggest that every observation and therefore every thickness estimate comes with the same level of uncertainty.

Here is another map from my thesis. It’s pretty dense, but shows the process in action:

The top panel is LIDAR ‘height above sea level’ for sea ice. Orange points are ‘sea level reference’ markers, and the grey patch highlights an intensive survey plot. The second panel is uncertainty associated with each height measurment. In panel  three we see modelled sea ice thickness (I’ll write another post about that later), and the final panel shows the uncertainty associated with each thickness estimate. Thickness uncertainties are greater than LIDAR height uncertainties because we’re also accounting for uncertainties in each of the other model parameters (LIDAR elevations are just one). So, when I get to publishing sea-ice thickness estimates, I can put really well made error bars around them!

ii. why you, as a commercial surveyor or an agency contracting a LiDAR survey or a LIDAR end user should bother:

The computation of uncertainties is straightforward and quick once the initial figuring of the Jacobians is done – and these only need to be recomputed when you change your instrument configuration. HeliMap (Switzerland) do it on-the-fly and compute a survey quality metric (see Schaer et al, 2007 [2]) which allows them to repeat any ‘out of tolerance’ areas before they land. Getting an aircraft in the air is hard, and keeping it there for an extra half an hour is easy – so this capability is really useful in terms of minimising costs to both contracting agencies and surveyors. This analysis in conjunction with my earlier post on ‘where exactly to measure in LiDAR heights‘ shows you where you can assign greater confidence to a set of LiDAR points.

It’s also a great way to answer some questions – for example are surveyors flying over GCP’s at nadir, and therefore failing to meet accuracy requirements off-nadir? (this is paranoid ‘all the world is evil’ me speaking – I’m sure surveyors are aware of this stuff and collect off-nadir ground control matches as well). Are critical parts of a survey being captured off-nadir, when it would be really useful to get the best possible data over them? (this has implications for flight planning). As a surveyor, this type of thinking will give you fewer ‘go back and repeat’ jobs – and as a contracting agency, you might spend a bit more on modified flight plans, but not a lot more to get actually really great data instead of signing off and then getting grief from end users.

As an end user of LiDAR products, If you’re looking for data quality thresholds – for example ‘points with noise < N m over a flat surface’, this type of analysis will help you out. I’ve also talked to a number of end users who wonder about noisy flight overlaps, and why some data don’t appear to be well-behaved. Again, having some quality metric around each point will help an end-user determine which data are useful, and which should be left alone.

Summarising

I certainly stand on the shoulders of giants here, and still incur brain-melting when I try to come at these concepts from first-principles (linear algebra is still hard for me!). However, the idea of being able to estimate in an a-priori quality metric is, in my mind, really useful.

I don’t have much contact with commercial vendors, so I can only say ‘this is a really great idea, do some maths and make your LiDAR life easer!’.

I implemented this work in MATLAB, and you can find it here:

https://github.com/adamsteer/LiDAR-georeference

With some good fortune it will be fully re-implemented in Python sometime this year. Feel free to clone the repository and go for it. Community efforts rock!

And again, refer to these works. I’ll see if I can find any more recent updates, and would really appreciate hearing about any recent work on this stuff:

[1] Glennie, C. (2007). Rigorous 3D error analysis of kinematic scanning LIDAR systems. Journal of Applied Geodesy, 1, 147–157. http://doi.org/10.1515/JAG.2007. (accessed 19 January 2017)

[2] Schaer, P., Skaloud, J., Landtwing, S., & Legat, K. (2007). Accuracy estimation for laser point cloud including scanning geometry. In Mobile Mapping Symposium 2007, Padova. (accessed 19 January 2017)

The LiDAR uncertainty budget I: georeferencing points

This is part 1 of 2, explaining how uncertainties in LiDAR point geolocation can be estimated for one type of scanning system. We know LiDAR observations of elevation/range are not exact (see this post), but a critical question of much interest to LiDAR users is ‘how exact are the measurements I have’?

As an end-used of LiDAR data I get a bunch of metadata that is provided by surveyors. One of the key things I look for are the accuracy estimates. Usually these come as some uncertainty in East, North and Up measurements, in metres, relative to the spatial reference system the point measurements are expressed in. What I don’t get is any information about how these figures are arrived at, or if they apply equally to every point. It’s a pretty crude measure.

As a LiDAR maker, I was concerned with the uncertainty of each single point – particularly height – because I use these data to feed a model for estimating sea ice thickness. I also need to feed in an uncertainty – so that I can put some boundaries around how good my sea ice thickness estimate is. However, there was no way of doing so in an off the shelf software package – so I implemented the LiDAR geoferencing equations and a variance-covariance propagation method for them in MATLAB  and used these. This was a choice of convenience at the time, and I’m now slowly porting my code to Python, so that you don’t need a license to make LiDAR points and figure out their geolocation uncertainties.

My work was based on two pieces of fundamental research: Craig Glennie’s work on rigorous propagation of uncertainties in 3D [1], and Phillip Schaer’s implementation of the same equations [2]. Assuming that we have a 2D scanner, the LiDAR georeferencing equation is:

\begin{bmatrix}  x \\  y \\  z \\  \end{bmatrix}^m=  \begin{bmatrix}  X\\  Y\\  Z\\  \end{bmatrix} + R^b_m  \left[  R^b_s  \rho  \left(  \begin{gathered}  sin\Theta\\  0\\  cos\Theta  \end{gathered}  \right)  +  \begin{bmatrix}  a_x\\  a_y\\  a_z\\  \end{bmatrix}^b  \right]

The first term on the right is the GPS position of the vehicle carrying the LiDAR:

\begin{bmatrix}  X\\  Y\\  Z\\  \end{bmatrix}

The next term is made up of a few things. Here we have points in LiDAR scanner coordinates:

\rho\left(  \begin{gathered}  sin\Theta\\  0\\  cos\Theta  \end{gathered}  \right)

…which means ‘range from scanner to target’ (\rho ) multiplied by sin\theta to give an X coordinate and cos\theta to give a Z coordinate of the point measured.

Note that there is no Y coordinate! This is a 2D scanner, observing an X axis (across track) and a Z axis (from ground to scanner). The Y coordinate is provided by the forward motion of a vehicle, in this case a helicopter.

For a 3D scanner, or a canner with an elliptical scan pattern, there will be additional terms describing where a point lies in the LiDAR frame. Whatever term is used at this point, the product is the position of a reflection-causing object in the LiDAR instrument coordinate system which is rotated to the coordinate system of the vehicle’s navigation device using the matrix:

R^b_s

The point observed also has a lever arm offset added (the distance in 3 axes between the navigation device’s reference point and the LiDAR’s reference point), so we pretend we’re putting our navigation device exactly on the LiDAR instrument reference point:

\begin{bmatrix}  a_x\\  a_y\\  a_z\\  \end{bmatrix}^b

This mess of terms is finally rotated to a mapping frame using euler angles in three axes (essentially heading, pitch, roll) recorded by the navigation device:

R^b_m

…and added to the GPS coordinates of the vehicle (which are really the GPS coordinates of the navigation system’s reference point).

There are bunch of terms there – 14 separate parameters which go into producing a LiDAR point, and that’s neglecting beam divergence, and only showing single returns. Sounds crazy – but the computation is actually pretty efficient.

Here’s a cute diagram of the scanner system I was using – made from a 3D laser scan and some engineering drawings. How’s that? Using a 3D scanner to measure a 2D scanner. Even better, the scan was done on the helipad of a ship in the East Antarctic pack ice zone!

You can see there the relationships I’ve described above. The red box is our navigation device – a dual-GPS, three-axis-IMU strapdown navigator, which provides us with the relationship between the aircraft body and the world. The green cylinder is the LiDAR, which provides us ranges and angles in its own coordinate system. The offset between them is the lever arm, and the orientation difference between the axes of the two instruments is the boresight matrix.

Now consider that each of those parameters from each of those instruments and the relationships between them have some uncertainty associated with them, which contributes to the overall uncertainty about the geolocation of a given liDAR point.

Mind warped yet? Mine too. We’re all exhausted from numbers now, so part 2 will examine how we take all of that stuff and determine, for every point, a geolocation uncertainty.

Feel free to ask questions, suggest corrections, or suggest better ways to clarify some of the points here.

There’s some code implementing this equation here: https://github.com/adamsteer/LiDAR-georeference – it’s Apache 2.0 licensed so feel free to fork the code and make pull requests to get it all working, robust and community-driven!

Meanwhile, read these excellent resources:

[1] Glennie, C. (2007). Rigorous 3D error analysis of kinematic scanning LIDAR systems. Journal of Applied Geodesy, 1, 147–157. http://doi.org/10.1515/JAG.2007. (accessed 19 January 2017)

[2] Schaer, P., Skaloud, J., Landtwing, S., & Legat, K. (2007). Accuracy estimation for laser point cloud including scanning geometry. In Mobile Mapping Symposium 2007, Padova. (accessed 19 January 2017)

LiDAR thoughts – where to measure, exactly?

LiDAR is a pretty common tool for geospatial stuff. It means ‘Light Detection and Ranging’. For the most part involved shining a laser beam at something then measuring how long a reflection takes to come back. Since we approximate the speed of light, we can use the round trip time to estimate the distance between the light source and ‘something’ with a great degree of accuracy. Modern instruments perform many other types of magic – building histograms of individual photons returned, comparing emitted and returned wave pulses, and even doing this with many different parts of the EM spectrum.

Take a google around about LiDAR basics – there are many resources which already exist to explain all this, for example https://coast.noaa.gov/digitalcoast/training/intro-lidar.html.

What I want to write about here is a characteristic of the returned point data. A decision that needs to be made using LiDAR is:

Where should I measure a surface?

…but what – wait? Why is this even a question? Isn’t LiDAR just accurate to some figure?

Sort of. A few years ago I faced a big question after finding that the LiDAR I was working on was pretty noisy. I made a model to show how noisy the LiDAR should be, and needed some data to verify the model. So we hung the LiDAR instrument in a lab and measured a concrete floor for a few hours.

Here’s a pretty old plot of what we saw:

What’s going on here? In the top panel, I’ve spread our scanlines along an artificial trajectory heading due north (something like N = np.arange(0,10,0.01)), with the Easting a vector of zeroes and height a vector of 3’s – and then made a swath map. I drew in lines showing where scan angle == 75, 90, and 115 are.

In the second panel (B), there’s a single scanline show across-track. This was kind of a surprise – although we should have expected it. What we see is that the range observation from the LiDAR is behaving as specified – accurate to about 0.02 m (from the instrument specifications). What we didn’t realise was that accuracy is angle-dependent, since moving away from instrument nadir the impact of angular measurement uncertainties becomes greater than the ranging uncertainty of the instrument.

Panels C and D show this clearly – near instrument nadir, ranging is very good! Near swath edges we approach the published instrument specification.

This left us with the question asked earlier:

When we want to figure out a height reference from this instrument, where do we look?

If we use the lowest points, we measure too low. Using the highest points, we measure too high. In the end I fitted a surface to the points I wanted to use for a height reference – like the fit line in panel B – and used that. Here is panel B again, with some annotations to help it make sense.

You can see straight away there are bears in these woods – what do we do with points which fall below this plane? Throw them away? Figure out some cunning way to use them?

In my case, for a number of reasons, I had to throw them away, since I levelled all my points using a fitted surface, and ended up with negative elevations in my dataset. Since I was driving an empirical model based on these points, negative input values are pretty much useless. This is pretty crude. A cleverer, more data-preserving method will hopefully reveal itself sometime!

I haven’t used many commercial LiDAR products, but one I did use a reasonable amount was TerraSolid. It worked on similar principles, using aggregates of points and fitted lines/planes to do things which required good accuracy  – like boresight misalignment correction.

No doubt instruments have improved since the one I worked with. However, it’s still important to know that a published accuracy for a LiDAR survey is kind of a mean (like point density)  – some points will have a greater chance of measuring something close to where it really is than others. And that points near instrument nadir are more likely to be less noisy and more accurate.

That’s it for now – another post on estimating accuracy for every single LiDAR point in a survey will come soon(ish).

The figure shown here comes from an internal technical report I wrote in 2012. Data collection was undertaken with the assistance of Dr Jan L. Lieser and Kym Newbery, at the Australian Antarctic Division. Please contact me if you would like to obtain a copy of the report – which pretty much explains this plot and a simulation I did to try and replicate/explain the LiDAR measurement noise.

ACE CRC, Airborne LiDAR and Antarctic sea ice

Between late June and late August 2015 I worked with the Antarctic Climate and Ecosystems Co-operative Research Centre (ACE CRC) to tidy up some long running loose ends with an airborne LiDAR project. This project is close to home – my PhD revolves around cracking some of the larger nuts associated with getting a science result from survey flights undertaken between 2007 and 2012. However, I’ve worked on one small subset of data – and the CRC was in need of a way to unlock and use the rest.

Many technical documents exist for the airborne LIDAR system, but the ‘glue’ to tie them together was lacking. In six weeks I provided exactly that. The project now has a strong set of documentation covering the evolution of the system, how navigate the myriad steps involved in turning raw logfiles from laser scanners, navigation instruments and GPS observations into meaningful data, and how to interpret the data that arise from the system. After providing a ‘priority list’ of flight data to work on, ACE CRC also took advantage of my experience to churn out post-processed GPS and combined GPS + inertial trajectories for those flights. The CRC also now has the tools to estimate point uncertainty and reprocess any flights from the ground up – should they wish to.

All of which means ACE CRC are in a position to make meaningful science from the current set of airborne LiDAR observations over East Antarctic sea ice.

Some of this – a part of my PhD work and a small part of the overall project – is shown here. A first-cut of sea ice thickness estimates using airborne LiDAR elevations, empirical models for snow depth, and a model for ice thickness based on the assumption that ice, snow and seawater all exist in hydrostatic equilibrium.

figure showing sea ice draft from LiDAR points
Sea ice draft estimated from airborne LiDAR elevations, an empirical model for snow, and a hydrostatic model for ice thickness from elevations and snow depth. This image shows a model of the underside of sea ice – blue and green points are elevations, purple to yellow points visible here are ice draft estimates. Why are some drafts ‘positive’? That’s a question currently being worked on…