Version 17 (modified by anch, 13 years ago) (diff)

--

Creating DEMs from ASTER or SRTM data

The ASTER GDEM is a global elevation dataset at 1-arc-second (30m) horizontal resolution, covering the earth from 83S-83N. More data about the dataset is available at http://www.ersdac.or.jp/GDEM/E/2.html, the dataset itself can be downloaded from http://asterweb.jpl.nasa.gov/gdem-wist.asp.

The Shuttle Radar Topography Mission mapped global elevations at a 3 arc second (90m posting) resolution between 60N and 56S. More information is available from http://srtm.usgs.gov/, and the data are available via the USGS Earth Explorer

Both of these are relatively coarse for ARSF DEM usage but are sufficient in some cases. To create a DEM from these data:

  1. Obtain DEM data. ASTER should be used in preference to SRTM unless there is a good reason not to.
  • SRTM from http://earthexplorer.usgs.gov. Expand the "Digital Elevation" link on the left and tick the box labelled "SRTM". Enter bounding co-ordinates for the area you want data for, clicking the green "+" icon after each one. Note that this currently seems not to always work correctly for locations in DMS format (at least in Firefox), you may need to switch to decimal degrees (which work fine). Click on the "Search" button in the bottom-right, and you should be presented with a list of items matching your search (there should be only one thing in it). Click on the "SRTM" link (you may need to enable pop-ups if you've got them blocked), and you will be presented with a list of matching SRTM tiles and download links - download the data you want to use. Note that only SRTM 3-arc-second resolution data is available for sites outside North America. For sites in North America higher resolution data can be obtained via the USGS Seamless Server
  1. Unzip the downloaded zip file. ASTER data will contain two GeoTIFFs (you want the "_dem" one) - the file name refers to the lat/lon of the bottom-left corner of the tile. For SRTM data, descend into the directory structure to find the file w001001.adf - this holds the actual data, though you can't move it without the rest of the directory structure.

If the ASTER data covers more than one tile, you will need to unzip lots of files - for convenience:

for zipfile in ls *.zip; do yes | unzip $zipfile; done

Each zipfile contains a file called Readme.pdf and piping yes into unzip lets us overwrite it each loop.

  1. Fire up Grass. Select a location in lat/long projection using WGS84 datum (create one if none available), make a note of the selected location and mapset names.
  1. Import the data file:

For ASTER data:

r.in.gdal input=<aster_tile>_dem.tif output=<MAP_NAME>

For SRTM data:

r.in.gdal input=w001001.adf output=<MAP_NAME>

For multiple tiles

If the data you have downloaded spans multiple tiles, then you will need to run the above but inside a for loop like so:

For ASTER data:

for tile in `ls *_dem.tif`;
do r.in.gdal input=$tile output=$tile.tile;
done

For SRTM data:

for tile in `ls *.adf`;
r.in.gdal input=$tile output=$tile.tile;
done

Once this is done, we need to stick all of the tifs together. First run

g.mremove rast=*

This will list all of the tiles you have just read in. The output should look something like:

Collecting map names for current mapset <Your_Mapset>...
The following files would be deleted:
g.remove rast=ASTGTM_N63W017_dem.tif.tile,ASTGTM_N63W018_dem.tif....

I will refer to the comma separated list after "rast=" as <TILELIST> from hereon - This is the bit you need.

Now run

g.region rast=<TILELIST>

And finally,

r.patch input=<TILELIST> output=<MAP_NAME>

Write down <MAP_NAME> somewhere - you will need it in a minute.

  1. Quit Grass and then start it up again. Select a location using UTM projection for the target area or create one if none is available - make sure you've got the right UTM zone. Note you can also use whatever other projection you want, but this guide assumes you want UTM.

Check http://www.gpsinformation.org/utm-zones.gif For the UTM zone

  1. Here we can use the manual method of calculating the eastings/northings or simply grab them from the tile(s) we have input.

Manual: Use a conversion utility (eg see the spreadsheet at http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM#Spreadsheet) to determine the boundaries of the target flight in meters within the selected UTM zone. If the flight spans two or more UTM zones, you will have to calculate the offsets appropriate to work out the east and west boundaries of the area. eg. If western boundary is zone 32 and eastern boundary is zone 33, for eastern value use zone 32, add 750000 to value from spreadsheet to find eastings relative to zone 32 bound.

  1. Set active region to calculated eastings/northings using
g.region n=<north val> s=<south val> e=<east val> w=<west val> res=50

For ASTER please put res=30

  1. Use r.proj to convert from lat/long to UTM:
    r.proj input=<MAP_NAME> location=<lat/long_location> mapset=<lat/long_mapset> output=<output_map_name>
    

Where <MAP_NAME> is what we created in step 4.

  1. If r.proj gives an error saying that the input map is outside the bounds of the current region, go back to step 6 and check your numbers

Other method: Set the northings/eastings to something very large e.g.

g.region n=50000000 s=0 e=50000000 w=-50000000 res=50

For ASTER please put res=30

  1. Use r.proj to convert from lat/long to UTM:
    r.proj input=<MAP_NAME> location=<lat/long_location> mapset=<lat/long_mapset> output=<output_map_name>
    

Where <MAP_NAME> is what we created in step 4.

9.

Set the new region

g.region rast=<output_map_name>

You can run

gis.m &

To run a GUI and check whether the DEM covers the correct area, or to make the area smaller.

10 Output ASCII DEM from converted map as normal (see Creation of DEMs from NextMap data)