Changes between Version 33 and Version 34 of Help/LeicaLidarDems


Ignore:
Timestamp:
Apr 19, 2017, 2:23:54 PM (8 years ago)
Author:
dac
Comment:

Removed all GRASS text and replaced with new DEM generation method.

Legend:

Unmodified
Added
Removed
Modified

  • Help/LeicaLidarDems

    v33 v34  
    11= Creating a DEM from a LIDAR point cloud =
    22
    3 '''If you are looking for an automated way to create a DEM from ARSF LiDAR data, in particular for use in APL, the [https://github.com/pmlrsg/arsf_dem_scripts ARSF DEM scripts] are recommended.'''
     3One of the most common uses of LiDAR data is to generate a high resolution Digital Elevation Model (DEM) raster. This page explains the steps used by NERC-ARF to generate DEMs from the point clouds provided as part of the standard delivery.
    44
    5 A tutorial using a previous version of the scripts is available [wiki:Help/DEM_scripts here]. However, unless you are alredy familiar with these using the new scripts is recommended.
     5== LiDAR only DEM ==
    66
    7 This page is designed more as a tutorial for making a DEM in GRASS.
    8 
    9 This page contains instructions for creating a simple DEM from the ASCII point clouds acquired from the Leica ALS50 II in 2009 onwards. Instructions are given for using your ARSF data in the GRASS GIS system.
    10 
    11 There is also a section giving instructions on how to make your DEM suitable for use in the azgcorr software.
    12 
    13 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]).
    14 
    15 ----------
    16 == Using GRASS ==
    17 The 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).
    18 
    19 When this is done the ASCII point data can be loaded in. 
    20 
    21 === Import the LiDAR ===
    22 
    23 To import the ASCII data it must be within the GRASS region limits, any data outside the current region will not be imported.
    24 
    25 '''Scan through the point cloud and find the min/max Eastings and Northings''':
     7For LiDAR deliveries from 2016 onwards the LiDAR DEM is generated with the NERC-ARF DEM scripts (https://github.com/pmlrsg/arsf_dem_scripts) using a command similar to:
    268
    279{{{
    28 r.in.xyz -s input=<LDRfilename> output=<outputmapname> x=2 y=3 z=4 fs=' '
     10las_to_dsm.py --method points2grid --projection UKBNG -r 2 \
     11    -o dem/RG15_18-2015_290a-lidar-bng.dem \
     12     flightlines/las1.2/*.LAS
    2913}}}
    3014
    31 where <LDRfilename> is the filename you want to read in, <outputmapname> is the name you wish to call this within GRASS, x, y, and z are equal to the column numbers which contain the Easting, Northing and elevation values, fs is the field separator. For more info on this command see the GRASS manual [http://grass.fbk.eu/gdp/html_grass63/r.in.xyz.html]
     15Where the projection is given in the readme file.
    3216
    33 '''Set the region such that it contains all the point cloud data, and the resolution you wish to use''' – in this case 2.0m:
     17This command will merge all LAS files together using [http://www.cs.unc.edu/~isenburg/lastools/download/lasmerge_README.txt `lasmerge`], dropping any points with class 7 (noise).
     18
     19To generate to DEM the [http://www.opentopography.org/otsoftware/points2grid points2grid] program is called as part of the script with the following options:
    3420
    3521{{{
    36 g.region  n=max_northing s=min_northing w=min_easting e=max_easting res=2.0
     22points2grid --exclude_class 7 --output_file_name /tmp/RG15_18-2015_290a-lidar-bng \
     23      --search_radius 2 --output_format arc --resolution 2 \
     24      --idw --fil_window_size 7 \
     25      --first_return_only
    3726}}}
    3827
    39 replacing the keywords  max_northing,  min_northing,  min_easting,  max_easting with the values from the `r.in.xyz` command above. You may wish to add a small buffer on to ensure points on the boundary are imported.
     28Points to grid uses inverse distance weighted interpolation to generate the DEM, more details about the method used are available from http://www.opentopography.org/otsoftware/points2grid. As only first returns are considered the resulting raster is a Digital Surface Model (DSM), representing the top of canopy and buildings.
    4029
    41 Now we can '''import the laser point cloud data into GRASS'''. This uses the same `r.in.xyz` command as above but without the -s flag:
     30To generate a LiDAR only DEM using the same method as deliveries prior to 2016 use the flag `--method GRASS`, this will take the average of all non-noise first return points within each grid cell.
     31
     32== Patched DEM for use in APL ==
     33
     34For processing hyperspectral data in APL a second DEM is supplied with delivered data, reprojected to WGS84 Lat/Long and patched with another DEM (normally ASTER) to extend coverage and fill in any gaps. For deliveries after 2016 a command similar to the one below is used  to generate a patched LiDAR and ASTER DEM:
    4235
    4336{{{
    44 r.in.xyz input=<LDRfilename> output=<outputmapname> x=2 y=3 z=4 fs=' '
     37create_dem_from_lidar.py --in_projection UKBNG --out_projection WGS84LL \
     38        --aster --gridded \
     39        -o dem/RG15_18-2015_290a-lidar_ASTER-wgs84_latlong.dem \
     40        RG15_18-2015_290a-lidar-bng.dem
    4541}}}
    4642
    47 The above 3 steps need to be repeated for each of the point cloud files you wish to use to make the DEM with.
     43For deliveries before 2016 the command was similar but a directory containing LAS files is passed in and the `--las` flag is used to indicate the LiDAR data should be gridded before patching.
    4844
    49 When all point cloud files are loaded into GRASS, the '''region needs to be changed such that it covers the area of all the point clouds''':
    50 
    51 {{{
    52 g.region rast=mapname1,mapname2,....,mapnameN
    53 }}}
    54 
    55 where mapname1, mapname2, etc are the map names of the imported point cloud data from the above steps.
    56 
    57 All the separate '''LIDAR maps can now be patched together to make a single large raster''' data set. To do this, use the following command:
    58 
    59 {{{
    60 r.patch in=mapname1,mapname2,... out=lidar_mosaic
    61 }}}
    62 
    63 where the map names are as before, the imported point cloud rasters, and lidar_mosaic is the output name for the single concatenated raster set. For more details on this command see [http://grass.fbk.eu/gdp/html_grass63/r.patch.html].
    64 
    65 
    66 === Dealing with small holes ===
    67 
    68 If '''small''' holes are present in your data, say due to small lakes or removed noise, the easiest way to fill these in is to patch the data with a lower resolution version of the LiDAR. This can be done by changing the region to a lower resolution and resampling:
    69 
    70 {{{
    71 g.region res=5
    72 r.resamp.stats input=lidar_mosaic output=lowres_mosaic
    73 }}}
    74 
    75 which 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.
    76 
    77 This 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).
    78 
    79 {{{
    80 r.patch in=lidar_mosaic,lowres_mosaic  out=high_low_patched
    81 }}}
    82 
    83 
    84 === Dealing with larger holes - using 3rd party data ===
    85 
    86 If 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.
    87 
    88 As this process is very similar to extending the region of the DEM please see how to do it [here].
    89 
    90 
    91 === Dealing with larger holes - interpolating existing data ===
    92 
    93 If you do not wish to use 3rd party data or none exists for your area then it is possible it interpolate and extrapolate.
    94 
    95 This 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.
    96 
    97 {{{
    98 r.fillnulls input=raster_to_fill output=filled_raster tension=40. smooth=0.1
    99 }}}
    100 
    101 which 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.
    102 
    103 
    104 '''Alternative Method'''
    105 
    106 If 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.
    107 
    108 First 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:
    109 
    110 {{{
    111 g.region res=20
    112 r.resamp.stats input=raster_name output=mask_name
    113 }}}
    114 
    115 To assign the mask and reset the region resolution (in this example to 2m) use:
    116 
    117 {{{
    118 g.region res=2
    119 r.mask input=mask_name
    120 }}}
    121 
    122 To interpolate the data we will use 'r.surf.idw' which requires the data to be scaled first to retain precision.
    123 
    124 This command will scale the data by 100:
    125 {{{
    126 r.mapcalculator amap=raster_to_scale formula=100*A outfile=rasterx100 help=-
    127 }}}
    128 
    129 Run the interpolation command:
    130 {{{
    131 r.surf.idw input=rasterx100 output=interpolatedx100
    132 }}}
    133 
    134 And rescale the data back
    135 {{{
    136 r.mapcalculator amap=interpolatedx100 formula=0.01*A outfile=interpolated help=-
    137 }}}
    138 
    139 this results in a map where the internal holes have been filled, but the area outside the swath coverage remains unchanged.
    140 
    141 
    142 === Filtering the data ===
    143 
    144 It may be wise to filter / smooth the data to remove localised spikes. This can be done using the following command:
    145 
    146 {{{
    147 r.neighbors input=raster_to_filter output=filtered_raster method=median size=3 --overwrite
    148 }}}
    149 
    150 which 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].
    151 
    152 
    153 === Extending the coverage using ASTER data ===
    154 
    155 If 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 
    157 === Writing out the DEM ===
    158 
    159 The DEM can be written out from GRASS in various formats (see [http://grass.fbk.eu/gdp/html_grass63/raster.html] functions starting 'r.out').
    160 
    161 To write the data out as an ASCII grid (for example to use in azgcorr) you could use:
    162 
    163 {{{
    164 r.out.ascii input=raster_name output=output_filename null=0
    165 }}}
    166 
    167 which 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].
    168 
    169 To write out the data as an ENVI BIL file (for example to use in APL) you could use:
    170 
    171 {{{
    172 r.out.gdal format=ENVI type=Float32 input=raster_name output=output_filename nodata=9999
    173 }}}
    174 
    175 which will write out raster 'raster_name' to a 32-bit floating point ENVI file named 'output_filename', using 9999 to represent null data.
    176 
    177 '''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.
    178 
    179 
    180 
    181 ------------
    182 == Making the DEM suitable for azgcorr ==#azgcorrheader
    183 To 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):
    184 
    185 `or c r xm ym xx yx gi`
    186 
    187 where:
    188 
    189 or = row order of the data [0 if South to North, 1 if North to South]
    190 
    191 c  = number of columns
    192 
    193 r  = number of rows
    194 
    195 xm = minimum easting
    196 
    197 ym = minimum northing
    198 
    199 xx = maximum easting
    200 
    201 yx = maximum northing
    202 
    203 gi = grid size (spacing)
    204 
    205 An example header, for a DEM of 2000 rows x 2000 columns at 5m resolution, might be: 1 2000 2000 400000 850000 410000 860000 5
    206 
     45For more details on this command, and usage outside the NERC-ARF-DAN systems see the [https://github.com/pmlrsg/arsf_dem_scripts/blob/master/doc/source/tutorial_lidar.md#create-a-lidar--aster-dsm-for-use-in-apl tutorial].