| 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=-`''' |
| | 76 | 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). |
| | 77 | |
| | 78 | {{{ |
| | 79 | r.patch in=lidar_mosaic,lowres_mosaic out=high_low_patched |
| | 80 | }}} |
| | 81 | |
| | 82 | |
| | 83 | === Dealing with larger holes - using 3rd party data === |
| | 84 | |
| | 85 | 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. |
| | 86 | |
| | 87 | As 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 | |
| | 92 | If you do not wish to use 3rd party data or none exists for your area then it is possible it interpolate and extrapolate. |
| | 93 | |
| | 94 | 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. |
| | 95 | |
| | 96 | {{{ |
| | 97 | r.fillnulls input=raster_to_fill output=filled_raster tension=40. smooth=0.1 |
| | 98 | }}} |
| | 99 | |
| | 100 | 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. |
| | 101 | |
| | 102 | |
| | 103 | '''Alternative Method''' |
| | 104 | |
| | 105 | 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. |
| | 106 | |
| | 107 | 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: |
| | 108 | |
| | 109 | {{{ |
| | 110 | g.region res=20 |
| | 111 | g.resamp.stats input=raster_name output=mask_name |
| | 112 | }}} |
| | 113 | |
| | 114 | To assign the mask and reset the region resolution (in this example to 2m) use: |
| | 115 | |
| | 116 | {{{ |
| | 117 | g.region res=2 |
| | 118 | r.mask input=mask_name |
| | 119 | }}} |
| | 120 | |
| | 121 | To interpolate the data we will use 'r.surf.idw' which requires the data to be scaled first to retain precision. |
| | 122 | |
| | 123 | This command will scale the data by 100: |
| | 124 | {{{ |
| | 125 | r.mapcalculator amap=raster_to_scale formula=100*A outfile=rasterx100 help=- |
| | 126 | }}} |
| | 127 | |
| | 128 | Run the interpolation command: |
| | 129 | {{{ |
| | 130 | r.surf.idw input=rasterx100 output=interpolatedx100 |
| | 131 | }}} |
| | 132 | |
| | 133 | And rescale the data back |
| | 134 | {{{ |
| | 135 | r.mapcalculator amap=interpolatedx100 formula=0.01*A outfile=interpolated help=- |
| | 136 | }}} |
| | 139 | |
| | 140 | |
| | 141 | === Filtering the data === |
| | 142 | |
| | 143 | It may be wise to filter / smooth the data to remove localised spikes. This can be done using the following command: |
| | 144 | |
| | 145 | {{{ |
| | 146 | r.neighbors input=raster_to_filter output=filtered_raster method=median size=3 --overwrite |
| | 147 | }}} |
| | 148 | |
| | 149 | 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]. |
| | 150 | |
| | 151 | |
| | 152 | === Extending the coverage using ASTER data === |
| | 153 | |
| | 154 | 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. |
| | 155 | |
| | 156 | === Writing out the DEM === |
| | 157 | |
| | 158 | 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'). |
| | 159 | |
| | 160 | To write the data out as an ASCII grid (for example to use in azgcorr) you could use: |
| | 161 | |
| | 162 | {{{ |
| | 163 | r.out.ascii input=raster_name output=output_filename null=0 |
| | 164 | }}} |
| | 165 | |
| | 166 | 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]. |
| | 167 | |
| | 168 | To write out the data as an ENVI BIL file (for example to use in APL) you could use: |
| | 169 | |
| | 170 | {{{ |
| | 171 | r.out.gdal format=ENVI type=Float32 input=raster_name output=output_filename nodata=9999 |
| | 172 | }}} |
| | 173 | |
| | 174 | which 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 | |