DEM scripts

This page details and explains the use of the DEM creation scripts:

lidardem.sh
asterdem.sh
lidar_aster_dem.sh

The scripts listed are written in bash.

These scripts are used to create a DEM from LiDAR and elevation model tiles. The scripts use grass GIS to read in the elevation tiles and perform patching and various other operations on the data to create a DEM.

GRASS GIS

GRASS is a GIS (Geograohic Information System) and details of it can be found at: ​http://grass.fbk.eu/. There is a detailed api and usage guide, it is available for free download at: ​http://grass.fbk.eu/download/index.php, as the scripts are written to run on linux, depending on your build of fedora run the following to install grass, the scripts use grass63 by default:

To get started download the template grass database: grass_db_template.zip​

• Run the command:
  • unzip grass_db_template.zip
• To check that the grass download works run: grass63 -text grass_db_template/UKBNG/PERMANENT/
  • If it did not work successfully check the grass version that you have and the path to the grass template database
  • If it did work successfully press Ctrl-D or type exit to exit

Creating a DEM from ASTER data using asterdem.sh script

ASTER tiles

ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) is a sensor on the Terra satellite, ASTER is partly being used to obtain elevation maps. The ASTER tiles are 1 degree x 1 degree lat/lon. The ASTER version 1 have an approximate resolution of 30 metres, version 2 however has an approximate resolution of 15 metres.
15 metre resolution (left), 30 metre resolution (right)


ASTER tiles are in wgs84 lat/lon projection.

Downloading Tiles

The tiles can be downloaded from either of:

• ​http://reverb.echo.nasa.gov/reverb/
  • Create an account or log in
  • Tiles are downloaded by putting on an ftp site, the time period for this depends on the number of tiles, the fewer the faster
  • The best method to select tiles is by a bounding box, drag the box over the desired region
  • Type ASTER into the "Search Terms"
  • Select "ASTER Global Digital Elevation Model V002" in the "Select Datasets"
  • Then click "Search for Granules"
  • Click accept on the Data quality summary
  • Add all items to the cart
  • Click view items in cart
  • Click Order at the bottom of the page
  • An order items list will appear, click set: select the data usage, and agree to the terms, make sure "Use these values for all applicable order items" is selected, then save
  • Click proceed
  • Click "Submit Order" at the bottom of the page
  • A confirmation e-mail should be sent when the tiles are ready for download from an ftp site
• ​http://www.gdem.aster.ersdac.or.jp/
  • Create an account or log in
  • A tile count of under 100 can be downloaded immediately, over 100 the download is done in batches
  • Click on Search in the dialog box on the left
  • A page with a world map on it should have been loaded
  • Tiles can be selected by clicking on "Select tiles directly" or "Select tiles by polygon"
  • Select the desired tiles and click next
  • Check that your selection is valid, then click next
  • Select the purpose of the tile usage
  • Then download the tiles.

The files that are used for the elevation model are called ASTGTM*_dem.tif

Choosing Tiles

When creating an ASTER DEM it is important to choose the correct tiles to use, as this could impact on the speed at which the script runs, as if tiles with large gaps between them are chosen when grass fills the nulls it has a lot of nulls to fill, so will take a very long time. I suggest that a square region of tiles should be chosen, as this means that the script should not need to fill a large number of nulls. If you require tiles in the region (50.23 - 52.56, -2.34 - 1.23) lat/lon:

To calculate the tiles round the numbers down to get (50 - 52, -3 - 1), so you need:

ASTGTM2_N50W003_dem.tif
ASTGTM2_N50W002_dem.tif
ASTGTM2_N50W001_dem.tif
ASTGTM2_N50E000_dem.tif
ASTGTM2_N50E001_dem.tif
ASTGTM2_N51W003_dem.tif
ASTGTM2_N51W002_dem.tif
ASTGTM2_N51W001_dem.tif
ASTGTM2_N51E000_dem.tif
ASTGTM2_N51E001_dem.tif
ASTGTM2_N52W003_dem.tif
ASTGTM2_N52W002_dem.tif
ASTGTM2_N52W001_dem.tif
ASTGTM2_N52E000_dem.tif
ASTGTM2_N52E001_dem.tif

The name of the tiles relates to the bottom left corner of the region so ASTGTM2_N50W003_dem.tif covers (50-50.9999, 3-3.0001).

Geoid Height Grid

A geoid spheroid height model is required to elevate the aster tiles to the correct geoid. As the tiles are in lat/lon it is best to elevate it to wgs84, because this is a good standard. To elevate the aster we use a geoid height grid available from:

• ​http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/ww15mgh.grd.z

When the file has been downloaded decompress it by running the command on the command line:

• gunzip ww15mgh.grd.z

Then to get the file in a format compatible with grass run, this gives a header compatible with grass:

• sed "s/ -90.000000 90.000000 .000000 360.000000 .250000 .250000/north: 90.000000 \nsouth: -90.000000 \nwest: 0.000000 \neast: 360.000000\nrows: 721\ncols: 1441/g" ww15mgh.grd -i

Separation Files

If projecting the aster data to another projection with a differing elevation model a separation file needs to be used to elevate the aster correctly. For example if projecting the DEM to a osgb36 UKBNG projection from wgs84 lat/lon the elevation model is different. To elevate the aster we need a description of the height differences between the geoid and the spheroid, a file describing the differences between the two can be found at:

• ​http://www.ordnancesurvey.co.uk/docs/gps/ostn02-osgm02-files.zip
• Once the file has been downloaded run the command:
  • mkdir OSTN02_OSGM02files; unzip OSTN02_OSGM02files.zip -d OSTN02_OSGM02files
• Now run:
  • grass63 -text /PATH/TO/grass_db_template/UKBNG/PERMANENT/
• Grass should have opened if it did not look at the section above on Grass GIS
• First get the region that is going to be read in:
  • r.in.xyz input=OSTN02_OSGM02_GB.txt output=OSTN02_OSGM02_GB.txt x=2 y=1 z=6 fs="," --overwrite -s
    • This should return:

Range:	min	max
x:	0.000000	700000.000000
y:	1.000000	876951.000000
z:	-81.603000	0.000000

• Then actually set the region based on this, include a buffer size of 1 metre, resolution is 1000 for this file:
  • g.region e=700001 w=-1 n=876952 s=0 res=1000
• Now read in the separation file:
  • r.in.xyz input=OSTN02_OSGM02_GB.txt output=OSTN02_OSGM02_GB.txt x=2 y=1 z=6 fs="," --overwrite
• Now to check that the reading in of the file has worked run:
  • gis.m
    • In the manager click on "Add raster layer"
    • Then click on "Base map" and select "OSTN02_OSGM02_GB.txt"
    • In the "Map Display" click on "Redraw all layers".
    • You should now see the raster map that has been read in
• Now the separation file has been read in a separation DEM can be created
  • The raster map can be output to a file with the following command:
    • r.out.gdal input=OSTN02_OSGM02_GB.txt output=uk_separation_file_UKBNG.dem format=ENVI type=Float32 nodata=0 --overwrite
    • rm -f uk_separation_file_UKBNG.dem.aux.xml
  • This DEM can be used as a separation file when projecting from osgb36 UKBNG to wgs84 lat/lon, however, we want it for projecting the other way. The DEM we have created here can be used later for use with the LiDAR, so do not remove it.
• To project the DEM:
  • First we need to get a vector representing the entire region of the DEM:
    • v.in.region output=vector_box type=area --overwrite
  • Then to get the resolution projected:
    • Save the region so that we can get back here:
      • g.region save=actualregion --overwrite
    • Set the region to the the bottom left corner based on the resolution:
      • g.region n=s+`echo \`g.region -p | grep nsres\` | awk '{print $2}'` e=w+`echo \`g.region -p | grep ewres\` | awk '{print $2}'`
    • Create a vector box based on the resolution:
      • v.in.region output=vector_box_res type=area --overwrite
    • Reset the region:
      • g.region actualregion
  • To get the DEM in another projection:
    • Project the resolution vector to the wgs84 lat/lon location:
      • v.proj input=vector_box_res location=UKBNG mapset=PERMANENT output=vector_box_res --overwrite
    • Set the region to be the size of the vector and then set the resolution:
      • g.region vect=vector_box_res rows=1 cols=1
    • Project the total region vector to the wgs84 lat/lon location:
      • v.proj input=vector_box location=UKBNG mapset=PERMANENT output=vector_box --overwrite
    • Set the total region:
      • g.region vect=vector_box
    • Project the separation file to the current location:
      • r.proj input=OSTN02_OSGM02_GB.txt location=UKBNG mapset=PERMANENT output=OSTN02_OSGM02_GB.txt --overwrite
    • Output the separation file:
      • r.out.gdal input=OSTN02_OSGM02_GB.txt output=uk_separation_file_WGS84LL.dem format=ENVI type=Float32 nodata=0 --overwrite
      • rm -f uk_separation_file_WGS84LL.dem.aux.xml

Configuring asterdem.sh

To run asterdem.sh on your system the files listed above need to be configured in the script to point to the correct places. For example, the template grass database in the script points to a place in the arsf directory and this needs to be changed to point to where your template grass database is located. Variables that need to be changed:

• TEMPLATE_GRASS_DATABASE=/users/rsg/arsf/usr/share/grass_db_template/
• GEOID_SPHEROID=/users/rsg/arsf/dems/geoid-spheroid/ww15mgh.grd
• WGS84LL_SEPARATION_FILE=/users/rsg/arsf/dems/aster/separation_files/uk_separation_file_WGS84LL.dem

Using asterdem.sh

asterdem.sh

Generates a binary lat-long ASTER DEM from downloaded raw ASTER data

Arguments:

-o: Final output filename for DEM
(Optional)

-h: Display usage and exit
-a: Directory containing ASTER tiles only. If not given then script

calculates tiles to use

-n: Navigation file to use, sbet or sol

Examples:

• Create an ASTER DEM from tiles in the directory aster_tiles and create an output file called a_dem_filename.dem:
  • asterdem.sh -a aster_tiles/ -o a_dem_filename.dem

Creating a DEM from LiDAR data using lidardem.sh script

Configuring lidardem.sh

To run lidardem.sh on your system the files listed above need to be configured in the script to point to the correct places. For example, the template grass database in the script points to a place in the arsf directory and this needs to be changed to point to where your template grass database is located. Variables that need to be changed:

• TEMPLATE_GRASS_DATABASE=/users/rsg/arsf/usr/share/grass_db_template/
• WGS84LL_SEPARATION_FILE=/users/rsg/arsf/dems/aster/separation_files/uk_separation_file_WGS84LL.dem
• GRASS_VERSION=grass63

I would advise that the grass version be 6.3 as the script was written to work with this grass version.

Creating a DEM from LiDAR and ASTER data using lidar_aster_dem.sh script

Configuring lidar_aster_dem.sh

To run lidar_aster_dem.sh on your system the files listed above need to be configured in the script to point to the correct places. For example, the template grass database in the script points to a place in the arsf directory and this needs to be changed to point to where your template grass database is located. Variables that need to be changed:

• TEMPLATE_GRASS_DATABASE=/users/rsg/arsf/usr/share/grass_db_template/
• WGS84LL_SEPARATION_FILE=/users/rsg/arsf/dems/aster/separation_files/uk_separation_file_WGS84LL.dem

Using lidar_aster_dem.sh

This is a wrapper script for asterdem.sh and lidardem.sh, so the options are similar to those above. The script is configured for internal use, so here is a out of house usage: lidar_aster_dem.sh
Generates a binary LiDAR ASTER DEM

Arguments:

-l: File or Directory containing LiDAR data
-z: Specify an input projection of the LiDAR files
-o: Final output filename for DEM, defulats to leaving in grass database
(Optional)

-h: Display usage and exit
-j: Create a jpg image of the DEM, followed by filename
-t: Whether to use ASCII files to create the DEM, defaults to using LAS

files. It is suggested to use this option if the las file are very
large

-b: Whether to calculate buffer, defaults to using NSEW buffer size of

2km (Projection may vary)

-p: Projection to output DEM in:

input proj4 string or UKBNG or ll or UTM[zero_padded_number](N|S)

(Power Users)

-c: Whether to use points of classification 7, defaults to just points

of classification 1

-a: Whether to use all returns, defaults to using first return only
-r: Resolution to input LiDAR files with
-R: Resolution of output DEM, defaults to 2
-g: Grass database to use, defaults to creating a GRASS database locally
-i: Navigation file to use
-T: Whether to look for 30 metre tiles aswell as 15
-n: North buffer, defaults to 2000
-s: South buffer, defaults to 2000
-e: East buffer, defaults to 2000
-w: West buffer, defaults to 2000
-A: Directory containing raw unzipped ASTER data. If not given then

script needs to be run within a project

-S geoid-spheroid separation file
-B: Comma separated North,South,East,West bounds to calculate buffer size

from

example:

lidar_aster_dem.sh -l las1.0/ -z UKBNG -o a_dem_filename.dem

Examples:

• Create a LiDAR ASTER DEM from tiles in directory aster_tiles and UKBNG las files in las1.0 and output to a file called dem_filename.dem (Uses default buffer size of 2000 in all directions
  • lidar_aster_dem.sh -l las1.0 -A aster_tiles -o a_dem_filename.dem
• Create a LiDAR ASTER DEM by calculating ASTER tiles and utm30N las files in las1.0 and output to a file called dem_filename.dem (Uses default buffer size of 2000 in all directions
  • lidar_aster_dem.sh -l las1.0 -z UTM30N -A aster_tiles -o a_dem_filename.dem
• Create a LiDAR ASTER DEM from tiles in directory aster_tiles and las file las_file.LAS and output to a file called dem_filename.dem, use buffer size of 100 in all directions
  • lidar_aster_dem.sh -l las_file.LAS -A aster_tiles -N 100 -S 100 -E 100 -W 100 -o a_dem_filename.dem