Changes between Version 28 and Version 29 of Help/LeicaLidarDems

Feb 28, 2012, 1:08:25 PM (12 years ago)



  • Help/LeicaLidarDems

    v28 v29  
    77'''This assumes that you have removed any noisy points that you do not wish to be used in the making of your DEM (see [wiki:Processing/LidarNoisyPoints]).'''
     9Please refer to this page on [wiki:Help/DEM_scripts using scripts for creation of DEMS] as it is a fully scripted approach. This page is designed more as a tutorial for making a DEM in GRASS.
    912== Using GRASS ==
    1013The first step is to create a location of the required area. See the GRASS reference manual for help with this. Alternatively you can download a template grass database from here: [raw-attachment:wiki:Help/ GRASS template database]. The template database contains both a UK National Grid and WGS84 Latitude/Longitude region. Start GRASS using a region suitable for the LiDAR data (creating one if it doesn't exist).
    12 When this is done the ASCII point data can be loaded in.  To import the ASCII data it must be within the GRASS region limits, any data outside the current region will not be imported.
     15When this is done the ASCII point data can be loaded in. 
     17=== Import the LiDAR ===
     19To import the ASCII data it must be within the GRASS region limits, any data outside the current region will not be imported.
    1421'''Scan through the point cloud and find the min/max Eastings and Northings''':
    55 [[Image(mosaic_with_holes.png,width=100,align=center)]]
    56 {{{
    57 #!html
    58 <p style="text-align: center">Figure: Mosaic of point cloud data gridded at 2m resolution. Holes (no data) are shown as white</p>
    59 }}}
     62=== Dealing with small holes ===
    6974which in this example will create a new raster called lowres_mosaic at 5m resolution. Note that this is not suitable for larger holes in the data and other methods should be used in this case such as patching external data or interpolating.
    72 6.
    73 '''`r.mapcalculator amap=lidar_mosaic formula=100*A outfile=lidar_mosaicx100 help=-`'''
    75 '''` input=lidar_mosaicx100 output=lidar_mosaic_idwx100`'''
    77 '''`r.mapcalculator amap=lidar_mosaic_idwx100 formula=0.01*A outfile=lidar_mosaic_idw help=-`'''
    79 where the output raster, lidar_mosaic_idw, has been interpolated using an inverse distance weighted formula. As well as filling in holes within the lidar swath, this will also interpolate over empty parts of the GRASS region, resulting in possibly unrealistic data values.
    81 [[Image(interp_full_region.png,width=100)]]
    82 {{{
    83 #!html
    84 <p style="text-align: center">Figure: Interpolated point cloud data over full GRASS region.</p>
    85 }}}
    88 This can be improved if you wish by using a mask when interpolating the data. To create a mask, the easiest way is to use your input data at a low resolution. This creates a raster which covers the LIDAR swath but contains no holes in the data (assuming the resolution is selected low enough).  To import the data at a lower resolution , repeat steps 1-3 for each ASCII point cloud setting the `res` variable to a suitable value, e.g. 50.0 and outputting to a new map.  Then repeat steps 4 and 5 to create a raster covering the combined LIDAR swath.
    90 Then to use the mask:
    92 7. '''`r.mask input=<maskmapname>`'''
    94 and perform the interpolation step (with precision maintaining commands before and after as above):
    96 8. Re-set the resolution to 2 for the LiDAR before any calculations.
    97 '''`g.region res=2`'''
    99 9.
    100 '''`r.mapcalculator amap=lidar_mosaic formula=100*A outfile=lidar_mosaicx100 help=-`'''
    102 '''` input=lidar_mosaicx100 output=lidar_mosaic_idwx100`'''
    104 '''`r.mapcalculator amap=lidar_mosaic_idwx100 formula=0.01*A outfile=lidar_mosaic_idw help=-`'''
     76This can then be patched onto the original higher resolution LiDAR data such that the high resolution takes priority (as it is listed first in the 'in' option list).
     79r.patch in=lidar_mosaic,lowres_mosaic  out=high_low_patched
     83=== Dealing with larger holes - using 3rd party data ===
     85If there are large areas of no data then it is recommended that 3rd party DEMs are used to fill these. ASTER and SRTM DEMs have near global coverage and are free to download. This assumes using the ASTER data.
     87As this process is very similar to extending the region of the DEM please see how to do it [here].
     90=== Dealing with larger holes - interpolating existing data ===
     92If you do not wish to use 3rd party data or none exists for your area then it is possible it interpolate and extrapolate.
     94This can be done using the following command. Note that this can take a long time to complete depending on the amount of null data that is present in the raster and the resolution and size of the raster.
     97r.fillnulls input=raster_to_fill output=filled_raster tension=40. smooth=0.1
     100which will fill in the null values from 'raster_to_fill' and store the new data in 'filled_raster'. The tension and smooth values are used to define the spline interpolation. Note that this method will spam the screen with a lot of "Warnings" which GRASS says can be ignored. See [] for more information.
     103'''Alternative Method'''
     105If you wished to use an alternative method of interpolation, and not extrapolate outside the edges of the data, then you could follow a procedure similar to this.
     107First create a mask that covers the region of the data. The easiest way is to use your input data at a low resolution. This creates a raster which covers the LIDAR swath but contains no holes in the data (assuming the resolution is selected low enough).  To create lower resolution data, say at 20m, use:
     110g.region res=20
     111g.resamp.stats input=raster_name output=mask_name
     114To assign the mask and reset the region resolution (in this example to 2m) use:
     117g.region res=2
     118r.mask input=mask_name
     121To interpolate the data we will use '' which requires the data to be scaled first to retain precision.
     123This command will scale the data by 100:
     125r.mapcalculator amap=raster_to_scale formula=100*A outfile=rasterx100 help=-
     128Run the interpolation command:
     129{{{ input=rasterx100 output=interpolatedx100
     133And rescale the data back
     135r.mapcalculator amap=interpolatedx100 formula=0.01*A outfile=interpolated help=-
    106138this results in a map where the internal holes have been filled, but the area outside the swath coverage remains unchanged.
     141=== Filtering the data ===
     143It may be wise to filter / smooth the data to remove localised spikes. This can be done using the following command:
     146r.neighbors input=raster_to_filter output=filtered_raster method=median size=3 --overwrite
     149which will apply a median filter of kernel size 3 to the raster_to_filter, resulting in filtered_raster. More information on the filters that can be applied can be found at [].
     152=== Extending the coverage using ASTER data ===
     154If you wish to use the DEM for processing your hyperspectral data then it may be that you need to extend the coverage such that it covers the hyperspectral data. This can be done by patching on ASTER data. Information on this, and how to access ASTER data, can be found [wiki:Help/DEM_scripts#CreatingaDEMfromASTERdatausingasterdem.shscript here], along with how to offset the ASTER data from geoidal heights to spheroidal heights.
     156=== Writing out the DEM ===
     158The DEM can be written out from GRASS in various formats (see [] functions starting 'r.out').
     160To write the data out as an ASCII grid (for example to use in azgcorr) you could use:
     163r.out.ascii input=raster_name output=output_filename null=0
     166which will write out the raster 'raster_name' to the file 'output_filename' and set any null values to 0. To use this DEM the header will need to be converted into a style azgcorr understands, such as [#azgcorrheader here].
     168To write out the data as an ENVI BIL file (for example to use in APL) you could use:
     171r.out.gdal format=ENVI type=Float32 input=raster_name output=output_filename nodata=9999
     174which will write out raster 'raster_name' to a 32-bit floating point ENVI file named 'output_filename', using 9999 to represent null data.
     176'''IMPORTANT NOTE.''' If you wish to use the DEM in APL it '''MUST''' be in WGS84 Latitude / Longitude. If you wish to do this and you are not working in a WGS84 Latitude/Longitude region then you will need to reproject the data into that region before outputting.
    133 == Making the DEM suitable for azgcorr ==
     206== Making the DEM suitable for azgcorr ==#azgcorrheader
    134207To make a suitable ASCII DEM for use in the azgcorr software, the header information of the ASCII DEM file must be in a certain format. The required format is to have a header of one line (the first line of the file) with the DEM data following. The format is (as given by the azgcorr help):