Creating Maps From LIDAR Data
This is the way I (Glenn Horrocks, that is) am currently generating the LIDAR maps we are using at the Navshield, SAREX and search operations within NSW. This page is an attempt at documenting the method so others can use it, and hopefully improve on it. I am continuously developing these processes and improving them so they might change over time. Also I do not claim this method is efficient, elegant or optimised…. it is good enough and for now that will do.
Step 1: Setup Software
First, your computer. Your computer will need to be powerful, preferably a workstation class machine, running Windows 64-bit and at least 16GB RAM. Generating these maps is very computer intensive and if you only have a cheap laptop you do not have a hope of getting anything done. Also, the images you will generate are going to be huge so 32 bit systems are not going to handle it either. For instance the 22km x 11km map I generated for the Snowy Mountains required about 30 hours processing time.
Software you will need:
- Karttapullautin from http://www.routegadget.net/karttapullautin/ (What the *&% does that mean? It is Finnish for automatic map generator). Get the 64 bit version.
- NSWtopo from https://github.com/mholling/nswtopo. This is an excellent package written by Matthew Hollingworth, who I understand is an ACT rogainer.
- NSWtopo requires several packages. Go through the webpage on this carefully, but you will need Ruby, GDAL, ImageMagick and Google Chrome. You will want to make KMZ maps so get 7-zip (https://www.7-zip.org/) as well. I do not use the other optional software (SQLite, pngquant, Inkscape).
- LAStools from https://rapidlasso.com/lastools/. This is required for pre-processing the LIDAR data for Karta.
I also use python as a scripting language to control things. Get the latest version of python from https://www.python.org/ (python 3.9 as of 2021) in a 64 bit version. You will need a few modules, imageio, numpy, matplotlib and simplekml. To install the numpy module the easiest way is to go to Christoph Gohlke’s python extension page at https://www.lfd.uci.edu/~gohlke/pythonlibs/ and download the numpy+mkl WHL file for your python version and OS. You install it with the command line “pip install <Path to numpy***.whl file>”. After numpy is installed you can install all the other modules with command line “pip install imageio matplotlib simplekml”. Note the command line will need to have administrator privilege (right click on START/Windows System/Command Line, and select More/Run as administrator).
You will also need to sort out the PATHs so your system can find the executables. I trust you can do this, you are not going to get far if this is beyond you.
Step 2: Downloading LIDAR Data
- Go to ELVIS (https://elevation.fsdf.org.au/). Make sure you are wearing your funky white jump suit.
- Zoom the map to cover the area you want to generate the map for.
- The ELVIS interface is great, and gives you access to a whole bunch of datasets. If you click on the “Layers” button on the green menu bar you can see the available options.
- Select the area you wish to download the data for your map. I use the “Draw a Bounding Box” option next to “Select area by:”.
- Once you have selected your area ELVIS will have a little think, and then come up with a window titled “Found XX datasets”.
- There will be some data from Geoscience Australia. I do not use this data.
- The data I use is titled NSW Government Spatial Services. It should show two main data sets:
- DEMs (Digital Elevation Model) are used by NSWtopo to generate contours. The 2 Metre or 5 Metre data is best, but you can use the 1 Metre data if you like (but the download will be much bigger for no benefit for most maps). The DEMs are in 2km x 2km tiles so you probably need to download a number of them. Note that you make things hard for yourself if you use both more than one DEM resolution in a map, NSWtopo cannot handle mixed resolution. So get a single resolution over your entire map if you can.
- Point Clouds are the data from the LIDAR system with only a minimal amount of processing. They are used by Karttapullautin to generate the maps. The files are in 2km x 2km tiles so you probably need a number of them. You will find LIDAR data from before about 2014 is of very poor quality, 2015 to 2018 is much better and 2018 to 2021 is of excellent quality, so if there are multiple datasets covering your area get the most recent one.
- Select all the datasets you want (both DEM and Point Clouds) and click the blue Download button. The download for a 10km x 10km is going to be several GB, so make sure your download allowance can handle it, you have enough disk space free and your internet connection is fast enough it can complete the download before the sun turns into a white dwarf.
- ELVIS will send you an email with a download link. It can take ELVIS several minutes to generate this email, especially for large downloads.
- My experience of the downloads have been good, the ELVIS system appears stable and reliable so I have not had problems with the large downloads.
Step 3:Setting Up Karttapullautin
You will need to unzip the downloaded file. The open up the Point Cloud data and unzip all the data files.
In your Karta working directory for this map I set up a few files to do the processing:
- Make a directory called “FullMapTiles”
- Make a second directory called “VegetationTiles”
- MakeXYZMap.py – Generates a batch file which converts the Point Cloud (LAS) files to XYZ format, runs Karta, then sorts out the output into the FullMapTiles and VegetationTiles directories.
- MergeMap.bat – runs Karta to merged the processed tiles into a larger map
- MakeGeotiffs.bat – runs GDAL to convert the PNG generated by Karta into an industry standard georeferenced GEOTIFF file.
- pullauta.ini – settings for Karta
You might need to make some changes to these scripts, such as:
MakeXYZMap.py – High and low quality LIDAR data requires different processing. See comments in the file for details.
MakeGeotiffs.bat – The bit “+zone=56” means this map is in the UTM zone 56. If you are in a different zone you will need to change this.
pullauta.ini – I have set this to generate 5m contour lines with no form lines. While you can theoretically control the way Karta works by adjusting this file I find it very hard to use and it is challenging to do better than the defaults.
Step 4: Run Karttapullautin
Run MakeXYZMap.py to generate the batch file MakeXYZMap.bat.
Run MakeXYZMap.bat to get Karta to process all the files. This will take hours and possibly days. If you feel adventurous you can split the script to run several Karta processing tasks simultaneously (be aware that each Karta task needs to be in its own directory with its own ini file. If multiple Karta tasks are in the one directory they will overwrite each other’s temp files. And I do not find the built-in multiprocessing option in Karta very useful so I do not recommend using it.
…… hours pass…..
Finished! Now run MergeMap.bat to merge the tiles into a large map. If you try to generate too large a map it will give you weird error messages, so keep maps to approximately 10km x 10km. Maps larger than that should be made into multiple tiles. Also choose what output you want from Karta:
- All output – contours, vegetation, rocks/cliffs in “FullMapTiles”; OR
- Just vegetation in “VegetationTiles”.
Step 5: Setup NSWtopo
Make sure you read the documentation on the NSWtopo webpage, it describes most of the available features in this package.
NSWtopo is a command line driven application, it has not GUI. This means it is ideal to run using a standard script.
The file doMap.bat is the batch file I have used for a recent exercise. Stepping through the contents of the file:
- Line 1: Initialise using the lat/long coordinates defined. These are the corner points of the map you want. The easiest way to get these points is by using Google Earth, then select Tools/Options – Show Lat/Long – and select “Decimal Degrees”. If you put points at your map corners you can then get the result in the correct format.
- Line 2: Add the topographic data. This is the normal contour lines, roads, trails, watercourses, buildings and everything. In fact if you stop processing at this point you will have produced a standard 1:25000 topographic map.
- Lines 3-5: Remove topographic layers which are not required. The contours are removed because we are going to use the contours from Karta.
- Line 6: Add the grid, in this case a 1000m UTM grid.
- Line 7: Add the LIDAR map as a single raster layer.
- Lines 9-14: Add extra tracks from GPX files.
- Line 16: Render the maps at 300dpi, and output in PNG, SVG (vector graphics), GEOTIFF, kmz and georeferenced PDF formats.
If the LIDAR gods are smiling at you, you now have a shiny new LIDAR map.
In BSAR we use Garmin Oregons as our GNSS units. To put the LIDAR map on these units you need to import them as a kmz file. While the process described above does produce a kmz file it is not compatible with the Garmin Oregons and will not work. To produce maps for the Garmin Oregon you need to produce a kmz which complies with these requirements: https://support.garmin.com/en-GB/?faq=FtEncUXbaE0xE04yZ7gTq5
It is nice of Garmin to summarise the requirements – but they missed out one vital requirement. The KML/KMZ can define overlay images (which is what is used here to georeference the map image) by two methods, LatLonQuad or LatLonBox (see https://developers.google.com/kml/documentation/kmlreference for details). The Garmin Oregon ONLY accepts LatLonBox georeferencing.
The python file MakeOregonMap.py takes care of this and produces a compliant file in KML format. Note you will need to update line 18 to point to your NSWtopo GEOTIFF file.
Once you have generated the KML file, load it in Google Earth and check it is correct. If everything is in order then save the KML object (File/Save/Save Place As) and choose KMZ as file type. This will convert the KML file from MakeOregonMap.py into a KMZ file compatible with the Garmin Oregons.
You should now be able to upload the KMZ file to the Oregon (put it in the custom map directory), restart the device and the map should be there. It is a good idea to test the map has rendered properly before using it in an operation. It is easy to stuff something up in this complex work flow.
If you find any bugs in this method let me know and I will try to fix them. If you know a better way to do something then I want to hear it! I am always updating these processes and hopefully improving them so I will update this page from time to time with my latest method.
Alternatives to Karttapullautin
Karta is the most limiting part of this process. It can only produce 5m and 2.5m contours, it does not handle poor LIDAR data well, it is not easily tuned, I have no access to the source code and the outputs are limited. Alternatives are:
- Write my own – I have already started this but it will take a long time to develop. This would be ideal as we could output the LIDAR data as proper vector layers and this will make it play much better with other GIS systems.
- OCAD looks like it is a much better system and overcomes all the issues with Karta, but it is commercial software and I am too cheap to get a license: https://www.ocad.com/en/
- Note that existing data sets for vegetation density (like the SPOT5 vegetation data) are not tuned to give vegetation conditions for a person walking on the ground. Most existing datasets on vegetation are things like canopy density which is not very useful.
Alternatives to NSWtopo
Any GIS system could replace NSWtopo. In fact a high-end GIS system will allow much more control over what is on the map and where. I use NSWtopo as it is free and scriptable so the process is largely automated. Note that any alternative system means you will have to find a way of getting it into the Garmin Oregons, as they have their own unique requirements.