The Shuttle Radar Topography Mission (SRTM) is a tremendous resource because of its high accuracy and its consistency as a “snapshot” of Earth’s surface, acquired over just 10 days during February 2000. It is also freely available for anyone! SRTM was a specially designed mission of the Endeavour space shuttle, with a specially built sensor antenna, and aimed to map the Earth’s land topography using radar interferometry. As a result, during the 10 days of the mission, the land surface topography between latitudes 60°N and 56°S, or about 80% of Earth’s surface, was mapped with a final cell size of 1 arc-second (about 30 m) and absolute elevation errors under 10 m. More details about the mission can be found in the paper by Farr and collaborators [1] (which is open access and can be downloaded here).

SRTM data have been publicly available relatively soon after the data was processed. Three versions exist: version 1, the original sensor raw data; version 2, a processed version with data voids (arising from weak radar returns or other anomalies caused by, for example, deep valleys); and version 3, the void-filled version, compiled with the help of ASTER DEM. During many years, though, only data over the United States was available as V3 at the maximum resolution of 1 arc-second, whereas data for the rest of the world was offered as V2, 3 arc-second, or about 90 m pixel size. Fortunately, a recent happy decision by the US government has finally made available all data at full resolution for the whole globe. The process is gradual and is expected to be finished by next year.

Luckily, for our Landsat 8 example, the Bay of Naples is already available with 1 arc-second as of November 2014. We proceed then to download the data and prepare them. We start by opening the Earth Explorer site (http://earthexplorer.usgs.gov/) and choose the area of interest, either with the path/row or, as shown below, specifying a coordinate area with the map in the menu to the left:

srtm_000

we then choose the SRTM dataset in the menu to the left, under the tab ‘Data Sets’ by navigating the menu as follows: DigitalElevation -> SRTM -> SRTM 1 Arc-Second Global:

srtm_001

Clicking in “Results” will show us the 1° x 1° tiles in which SRTM data are compiled. We simply download them, using the icon with a green arrow next to each result item in the results list. I here chose to get the four westernmost tiles:

srtm_002

Once downloaded we get the four files:


n40_e013_1arc_v3.tif
n40_e014_1arc_v3.tif
n41_e013_1arc_v3.tif
n41_e014_1arc_v3.tif

These are GeoTIFF files. It is also possible to get the raw format grids, “hgt” files, but these require a bit more involved processing, so we here proceed with the GeoTIFFs. The name of the files refer to the coordinates of the lower left corner of each tile. Each of the files contains 3601 x 3601 elements, in 16 bit integer data type, with missing values signaled as -32767. The file coordinates are given as longitude and latitude, referred to the WGS84 ellipsoid, with elevation expressed in meters referred to the mean sea level (EGM96 geoid). As is tradition in this blog, we use the free and open source tools GDAL and Generic Mapping Tools (GMT) from the command line:

> gdalinfo n40_e013_1arc_v3.tif
Driver: GTiff/GeoTIFF
Files: n40_e013_1arc_v3.tif
Size is 3601, 3601
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (12.999861111111111,41.000138888888884)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
AREA_OR_POINT=Point
DTED_CompilationDate=0002
DTED_DataEdition=02
DTED_DigitizingSystem=SRTM
DTED_HorizontalAccuracy=0009
DTED_HorizontalDatum=WGS84
DTED_MaintenanceDate=0000
DTED_MaintenanceDescription=0000
DTED_MatchMergeDate=0000
DTED_MatchMergeVersion=A
DTED_NimaDesignator=DTED2
DTED_OriginLatitude=0400000N
DTED_OriginLongitude=0130000E
DTED_Producer=USCNIMA
DTED_RelHorizontalAccuracy=NA
DTED_RelVerticalAccuracy=0011
DTED_SecurityCode_DSI=U
DTED_SecurityCode_UHL=U
DTED_UniqueRef_DSI=F02 011
DTED_UniqueRef_UHL=F02 011
DTED_VerticalAccuracy_ACC=0008
DTED_VerticalAccuracy_UHL=0008
DTED_VerticalDatum=E96
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 12.9998611, 41.0001389) ( 12d59'59.50"E, 41d 0' 0.50"N)
Lower Left ( 12.9998611, 39.9998611) ( 12d59'59.50"E, 39d59'59.50"N)
Upper Right ( 14.0001389, 41.0001389) ( 14d 0' 0.50"E, 41d 0' 0.50"N)
Lower Right ( 14.0001389, 39.9998611) ( 14d 0' 0.50"E, 39d59'59.50"N)
Center ( 13.5000000, 40.5000000) ( 13d30' 0.00"E, 40d30' 0.00"N)
Band 1 Block=3601x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Unit Type: m

We now proceed to compile a mosaic with our four files:

> gdal_merge.py -o bayofnaples.tif -n -32767 *.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

Note that we have specified the no-data value, with the option ‘-n’. We now have a file with the following characteristics:

> gdalinfo bayofnaples.tif
Driver: GTiff/GeoTIFF
Files: bayofnaples.tif
Size is 7201, 7201
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (12.999861111111111,42.000138888888884)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 12.9998611, 42.0001389) ( 12d59'59.50"E, 42d 0' 0.50"N)
Lower Left ( 12.9998611, 39.9998611) ( 12d59'59.50"E, 39d59'59.50"N)
Upper Right ( 15.0001389, 42.0001389) ( 15d 0' 0.50"E, 42d 0' 0.50"N)
Lower Right ( 15.0001389, 39.9998611) ( 15d 0' 0.50"E, 39d59'59.50"N)
Center ( 14.0000000, 41.0000000) ( 14d 0' 0.00"E, 41d 0' 0.00"N)
Band 1 Block=7201x1 Type=Int16, ColorInterp=Gray

Since we plan to use this elevation model with our image from the Vesuvius volcano (this one), we proceed to clip the DEM to the same area as the image, reprojecting it at the same time to the UTM33N projection, and taking care to specify a 30 m pixel size, so there will be a one-to-one correspondence between the cells in the image and DEM. We accomplish this using gdalwarp:

> gdalwarp -t_srs 'EPSG:32633' -te 430995 4489005 470985 4539015 -tr 30 30 bayofnaples.tif bayofnaples_UTM33N.tif
Creating output file that is 1333P x 1667L.
Processing input file bayofnaples.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.

our new file has the following properties:

> gdalinfo bayofnaples_UTM33N.tif
Driver: GTiff/GeoTIFF
Files: bayofnaples_UTM33N.tif
Size is 1333, 1667
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32633"]]
Origin = (430995.000000000000000,4539015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 430995.000, 4539015.000) ( 14d10'46.25"E, 40d59'57.85"N)
Lower Left ( 430995.000, 4489005.000) ( 14d11' 6.14"E, 40d32'56.15"N)
Upper Right ( 470985.000, 4539015.000) ( 14d39'17.97"E, 41d 0' 6.51"N)
Lower Right ( 470985.000, 4489005.000) ( 14d39'26.33"E, 40d33' 4.67"N)
Center ( 450990.000, 4514010.000) ( 14d25' 9.19"E, 40d46'32.18"N)
Band 1 Block=1333x3 Type=Int16, ColorInterp=Gray

Finally, we can display a nice color plot using GMT. We begin by translating our GeoTIFF to a NetCDF grid:

> gdal_translate -of NetCDF bayofnaples_UTM33N.tif bayofnaples_UTM33N.grd
Input file size is 1333, 1667
0...10...20...30...40...50...60...70...80...90...100 - done.

We then clip the values below 0, setting them to zero, to facilitate colour coding:

> grdclip bayofnaples_UTM33N.grd -Gnew.grd -Sb0/0

The new grid, new.grd, has the following properties:

> grdinfo new.grd
new.grd: Title: GDAL Band Number 1
new.grd: Command: grdclip bayofnaples_UTM33N.grd -Gnew.grd -Sb0/0
new.grd: Remark:
new.grd: Gridline node registration used [Cartesian grid]
new.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
new.grd: x_min: 431010 x_max: 470970 x_inc: 30 name: x coordinate of projection [m] nx: 1333
new.grd: y_min: 4489020 y_max: 4539000 y_inc: 30 name: y coordinate of projection [m] ny: 1667
new.grd: z_min: 0 z_max: 1442 name: GDAL Band Number 1
new.grd: scale_factor: 1 add_offset: 0
new.grd: format: netCDF-4 chunk_size: 134,129 shuffle: on deflation_level: 3

We then proceed to make a nice color table, based on the master table dem2, using makecpt:

> makecpt -Cdem2 -T0/1500/10 > color.cpt

For visualization, it might be nicer if we indicate that the lowest values, the sea level, have to be blue. We do this by opening a text editor and adding a first line that encompasses the values from 0 to 0.1, with blue color (0/0/255 stands for 0 red, 0 green, 255 blue):

0 0/0/255 0.1 0/0/255
0.1 5.2267/105.17/63.16 10 5.2267/105.17/63.16
10 15.68/121.5/47.48 20 15.68/121.5/47.48

We can then plot the grid using grdimage:

> grdimage new.grd -JX15 -Rnew.grd -Ccolor.cpt -X1 -Y1 -P > plot.ps

and export the postscript file as a raster, and clipping the original page to the area of the map:

> ps2raster -A plot.ps

Digital elevation model covering the same area that the Vesuvius image, and with the same pixel size.
Digital elevation model covering the same area that the Vesuvius image, and with the same pixel size.

Presto! In the next post we will compute grids of aspect and slope from this DEM.

[1] Farr, T.G. et al. 2007. The shuttle radar topography mission. Reviews of geophysics 45.2.

Advertisements