Changes between Version 28 and Version 29 of Help/LeicaLidarDems


Ignore:
Timestamp:
Feb 28, 2012 1:08:25 PM (7 years ago)
Author:
mark1
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • Help/LeicaLidarDems

    v28 v29  
    66
    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]).'''
     8
     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.
     10
    811----------
    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/DEM_scripts:grass_db_template.zip 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).
    1114
    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. 
     16
     17=== Import the LiDAR ===
     18
     19To import the ASCII data it must be within the GRASS region limits, any data outside the current region will not be imported.
    1320
    1421'''Scan through the point cloud and find the min/max Eastings and Northings''':
     
    5360
    5461
    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 ===
     63
     64[[Image(mosaic_with_holes.png,width=150,align=right,margin=20)]]
    6065
    6166
     
    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.
    7075
    71 
    72 6.
    73 '''`r.mapcalculator amap=lidar_mosaic formula=100*A outfile=lidar_mosaicx100 help=-`'''
    74 
    75 '''`r.surf.idw input=lidar_mosaicx100 output=lidar_mosaic_idwx100`'''
    76 
    77 '''`r.mapcalculator amap=lidar_mosaic_idwx100 formula=0.01*A outfile=lidar_mosaic_idw help=-`'''
    78 
    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.
    80 
    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 }}}
    86 
    87 
    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.
    89 
    90 Then to use the mask:
    91 
    92 7. '''`r.mask input=<maskmapname>`'''
    93 
    94 and perform the interpolation step (with precision maintaining commands before and after as above):
    95 
    96 8. Re-set the resolution to 2 for the LiDAR before any calculations.
    97 '''`g.region res=2`'''
    98 
    99 9.
    100 '''`r.mapcalculator amap=lidar_mosaic formula=100*A outfile=lidar_mosaicx100 help=-`'''
    101 
    102 '''`r.surf.idw input=lidar_mosaicx100 output=lidar_mosaic_idwx100`'''
    103 
    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).
     77
     78{{{
     79r.patch in=lidar_mosaic,lowres_mosaic  out=high_low_patched
     80}}}
     81
     82
     83=== Dealing with larger holes - using 3rd party data ===
     84
     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.
     86
     87As this process is very similar to extending the region of the DEM please see how to do it [here].
     88
     89
     90=== Dealing with larger holes - interpolating existing data ===
     91
     92If you do not wish to use 3rd party data or none exists for your area then it is possible it interpolate and extrapolate.
     93
     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.
     95
     96{{{
     97r.fillnulls input=raster_to_fill output=filled_raster tension=40. smooth=0.1
     98}}}
     99
     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 [http://grass.fbk.eu/gdp/html_grass63/r.fillnulls.html] for more information.
     101
     102
     103'''Alternative Method'''
     104
     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.
     106
     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:
     108
     109{{{
     110g.region res=20
     111g.resamp.stats input=raster_name output=mask_name
     112}}}
     113
     114To assign the mask and reset the region resolution (in this example to 2m) use:
     115
     116{{{
     117g.region res=2
     118r.mask input=mask_name
     119}}}
     120
     121To interpolate the data we will use 'r.surf.idw' which requires the data to be scaled first to retain precision.
     122
     123This command will scale the data by 100:
     124{{{
     125r.mapcalculator amap=raster_to_scale formula=100*A outfile=rasterx100 help=-
     126}}}
     127
     128Run the interpolation command:
     129{{{
     130r.surf.idw input=rasterx100 output=interpolatedx100
     131}}}
     132
     133And rescale the data back
     134{{{
     135r.mapcalculator amap=interpolatedx100 formula=0.01*A outfile=interpolated help=-
     136}}}
    105137
    106138this results in a map where the internal holes have been filled, but the area outside the swath coverage remains unchanged.
     139
     140
     141=== Filtering the data ===
     142
     143It may be wise to filter / smooth the data to remove localised spikes. This can be done using the following command:
     144
     145{{{
     146r.neighbors input=raster_to_filter output=filtered_raster method=median size=3 --overwrite
     147}}}
     148
     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 [http://grass.fbk.eu/gdp/html_grass63/r.neighbors.html].
     150
     151
     152=== Extending the coverage using ASTER data ===
     153
     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.
     155
     156=== Writing out the DEM ===
     157
     158The DEM can be written out from GRASS in various formats (see [http://grass.fbk.eu/gdp/html_grass63/raster.html] functions starting 'r.out').
     159
     160To write the data out as an ASCII grid (for example to use in azgcorr) you could use:
     161
     162{{{
     163r.out.ascii input=raster_name output=output_filename null=0
     164}}}
     165
     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].
     167
     168To write out the data as an ENVI BIL file (for example to use in APL) you could use:
     169
     170{{{
     171r.out.gdal format=ENVI type=Float32 input=raster_name output=output_filename nodata=9999
     172}}}
     173
     174which will write out raster 'raster_name' to a 32-bit floating point ENVI file named 'output_filename', using 9999 to represent null data.
     175
     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.
     177
     178
     179
    107180
    108181[[Image(interp_with_mask.png,width=100)]]
     
    131204
    132205------------
    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):
    135208