I have downlowded the complete set (872 files) of CGIAR / SRTM 90 meter data tiles and most if not all (2780) available ASTER GDEM 30 meter files onto Jupiter. After creating a set of small (~ 12 - 16 map tiles) terrain map mosaics at 90 meter and 1 kilometer spatial resolutions, I have launched the first 'production scale' script that will create a terrain map mosaic of the continental United States at 1 kilometer resoution, comprised entirely of CGIAR SRTM data.

(continuing from previous update)
Digital Terrain File data sources:

CGIAR / SRTM 90 meter Digital Elevation Database Ver 4.1 obtained
from the 'bulk download' and interactive links on this page:
http://www.cgiar-csi.org/data/elevation/item/45-srtm-90m-digital-elevation-database-v41

The 'bulk' download ZIP file contained 872 image tiles,in separate ZIP files; each ZIP file containing one tile in ESRI .ASC (ASCII) format. All the ZIP files require a total of 14.1 GB disk space.

I converted the ASC images to GeoTiff format using the command:

gdal_translate -ot Int16 -of GTiff  $iName $iNewName

The 872 GeoTIFF images require 60 GB disk space, which is approximately one-third the amount required by the unZIPPd .ASC files.

Note: At least one of the bulk download files was defective (srtm_25_13.zip), and I replaced it by ordering from the GoogleMaps interface found on the above link. Note that individual tile downloads are available as GeoTIFF OR ASC format.

ASTER Global Digital Elevation Model (GDEM): from NASA WIST Portal:

https://wist.echo.nasa.gov/api/ (account required to download data)

I intended to download complete ASTER GDEM coverage for the Northern Hemisphere North of 60 degrees Latitude. To do so, I submitted approximately eight 'Lat/Long' bounding box search requests, which together covered the Northern Hemisphere. Over the next 24 hours I received emails with links to the search results on the WIST server.

I downloaded 2779 .ZIP format image files; the zip files consume 155 GB of disk storage. Each Zip file contains two GeoTIFF format images: the Terrain Elevation file (.dem) and the Quality Assessment file (.num) - and an XML format statistics summary (.aux). The .dem and .num files are of equal size, and each (2779 member) set consumes 70 GB of disk space.

Creating fused terrain (mosaic) images:

After acquiring the required input data, I have begun to assemble fused mosaics containing multiple imput data tiles. I have created Unix shell scripts based
upon the GDAL utility program 'gdalwarp'. This is a powerful utility with many
command options that specify input and output file characteristics. Refer to:
http://www.gdal.org/gdal_utilities.html for detailed command summaries of gdalwarp,
gdal_translate, and other utilites.

The gdalwarp command line used to create initial image mosaics

gdalwarp -of GTiff -ot Int16 -r near  -tr .0090436 .0090436 -srcnodata -9999 -dstnodata -9999 \ <list of input and output files.....>

gdalwarp -of GTiff -ot Float32 -r bilinear -tr .0090436 .0090436 -srcnodata -9999 -dstnodata -9999 \ <list of input and output files.....>

Note that since both SRTM and ASTER data sets use geographic coordinate space, there is no need to specify projection parameters for either data set.

Note on output image pixel size: The ASTER GDEM data has a pixel size of 30 meters
(or 0.000277777777778 degrees); the SRTM GDEM data has a pixel size of 90 meters
(or 0.000833333333333 degrees). I am creating the initial output mosaics with one kilometer pixels, using a pixel size of .009043 degrees, based on the calculations
at the site:
http://www.digitalmobilemap.com/index.php?slug=latitude-and-longitude-getting-started -
which has a good discussion of the topic.

Using these calculations, I obtain slightly different pixel sizes for 30 and 90 meters:

30 meters: .00027131
90 meters: .00081393

....different from the values in the incoming data.
We need to select a 'standard' pixel size for the 90 meter and 1 kilometer output data layers.

To date, I have created two mosaics: A small test area covering the interface between the ASTER and SRTM data in central Canada (90 meter resolution), and an 'SRTM-only' mosaic covering the Western Hemisphere (North and South America) at 1 kilometer resolution. Next will be a more extensive ASTER/SRTM mosaic, at 1 kilometer resolution.

Thanks Rick. Regarding pixel sizes (cell resolution), what we've been calling "1km" is just shorthand for 30 arcseconds. So for the unprojected (geographic) data layers that will ultimately be produced under this project, this means cell resolutions of 30/(60*60), or 0.0083333333333. Similarly the "30m" ASTER product is really 1 arcsecond (0.000277777... degrees), and the "90m" SRTM product is really 3 arcsecond (0.000833333... degrees).

That said, I don't think you should be generating resampled 1km rasters (fused or otherwise) yet, particular as I thought we were leaning toward producing the derived terrain layers (slope, aspect, etc) at 90m (3 arcseconds), and then aggregating. This means the DEM fusion has to happen at 90m, not 1km, and moreover validation and further post-processing should happen at the finer resolution too -- unless, I suppose, you can definitively show that artifacts and e.g. boundary anomalies evident at that resolution "disappear" when scaling up, not only for the DEM itself but also for all derived terrain layers.

Also, when you're simply combining tiles, you can get big efficiency gains by using the VRT format with GDAL, rather than using gdalwarp to combine the separate tif tiles into a big merged tif. For example, this builds a VRT file that "combines" all SRTM tiles in the working directory (assuming consistent naming):

gdalbuildvrt srtm_merged.vrt srtm_*.tif

This can yield two big benefits. First, this VRT file really just stores info about the input tiles and how they should be combined. It's a "lazy" format, which means it doesn't actually do the merging until it needs to (i.e., until you try to use it at some later point). So this way we don't have to store a large amount of what is essentially duplicated data (unmerged tiles and merged raster). Then you can still use other GDAL commands on this vrt file as though it were the real combined raster, including e.g. producing a GeoTIFF that really does contain all the data:

gdal_translate -of GTiff -ot Int16 srtm_merged.vrt srtm_merged.tif

The second benefit is that, in my experience, this can be considerably faster than using gdalwarp.

I see your point re: cell resolution, if we stay in the realm of arc-seconds.
Yes, we are in agreement about producing the derived layers from 90 m data and then aggregating.
We wouldn't want to use the fused 1 KM DEM layer for the derived layer calculations.
And the VRT format will be useful for creating the derived layers.

I have moved the "ReevesSpatialSoftwareStack.txt" file to the "Documents" section of Redmine, where it will be more obvious and accessible than buried within this Issue. Please make future edits/additions to the file under the new document structure in the Documents section. Rick-- could you remove it from here so that we don't get mixed up about which is the current version. Thanks!

Here, for the team's comment and feedback, is my proposed plan for producing the global Fused 1 Km Terrain layer.

******************************************************************************

Comprehensive Plan: Develop Global Fused 1 Km Terrain Layer
Delivery Date: June 3, 2011
Step one: Acquire incoming DEM Data Files
   CGIAR / SRTM
   ASTER GDEM
   CANADA DEM (boundary edge issue resolution and verification)

   Inventory of files: Deterimine completeness of global coverage / determine 'no coverage' areas
      comprehensive file list / GIS display

   Status: 95% complete for Western Hemisphere

Step Two: Initial Data Review and exploration
   Objectives:
      Gain experience with all DEM data types used in project AND software tools
      Develop, test, refine batch jobs that use GDAL tools for processing large
      numbers of DEM files into image mosaics

Focus: British Columbia / 60 Degree North Latitude mosaic boundary
   Observe and measure magnitude of boundary artifact
   Evaluate Quality Assurance validation methods: Example: ASTER GDEM '.num' files
   Report results to project scientists

   Status: 95% complete

Step Three: Production planning
   Objective: Organize large volume of DEM mosaic development task into managable, monitor reproducable segments

   Subdivid global study area into four quadrants
   Proposed (based on experience in Steps One and Two:
   Western Hemisphere North and South
       North America pilot study area: develop/evaluate/test Derived Data Product production
   Eastern Hemisphere North and South

   Create Image Mosaic production  Bacth Job that uses GDAL commands to generate DEM mosaic images
   Includes 'standard set' of GDAL / gdalwarp command parameters used for production of all four quadrants
      Command parameters will be adapted in response to lessons learned, production changes.

   Create one batch job per production quadrant + North America 'pilot study' area - including appropriate DEM tile file names.

   Status: 10% complete

Step Four: Production Pilot using North America 'pilot region'.
   Objective: Create fused DEM terrain layer AND subset of Derived Data Products
      Proposed: Slope, Aspect, Terrain Wetness Index ***(Suggestions from scientists desired)
      Develop derived data product production scripts
      Identify and resolve fused artifacts in DEM terrain layer that affect quality of derived data products
      Evaluate Disk Storage and CPU requirement issues in anticipation of Production (Step Five)
   Generate North America DEM 1 KM terrain layer
   Generate 90 meter (3 arcsecond) terrain layer
   Generate selected derived data products
       *** Work with Jim R, Mark, and scientific team members to refine methods for generating derived products
            ***** Baseline/Starting Point: AML Scripts and methods crated by Dr. Gallant and Ming

   Develop, test, document Quality Assurance validation methods
      Candidate 'collateral data': ASTER GDEM .num files

   Continue to monitor and record production time parameters
   Distribute pilot study data files to project scientests for evaluation
   Incorporate their comments and results into adapted production process

   Status: Not begun

Step Five: Production of complete Fused 1 Km Terrain Layer
   Objective: Produce the complete global data set by promised deadline

   In sequence, execute the production batch jobs, monitor and evaluate results, adjust process as necessary
   Provide examples to project scientists for evaluation and testing

   Status: Not begun

Notice: The 'implementation details' of this plan will be encapsulated into the directory file structure,
file lists for the global quadrants, and batch job commands.

****************************************************************************** 
Looking forward to timely discussion and implementation of this plan 

--Rick 

Here is an example of the workflows that generate fused digital terrain data layers
for one of the four Earth 'quadrants' (Northwest,Southwest,Northeast,Southeast)

Workflow consists of three shell scripts. Each script consists of one 'gdalwarp' commnad
and a collection of image files (.tif or .IMG). Here are the scripts, residing on /Jupiter:

/data/project/organisms/rcr/OutProducts/WestHemi/NorthWest/

1) mosaicNwHemiAster30ArcSecBL.sh  - generates mosaic of ASTER GDEM tiles, Northwest Quadrant
2) mosaicNwHemiCgiar30ArcSecBL.sh  -        "            CGIAR SRTM   "          " 
3) mosaicNwHemiAsterCgiar30ArcSecBL.sh - makes Northwest quadrant mosaic from 1), 2) results 

The core of the script is the gdalwarp command. Here is an example from 3) 

gdalwarp -of HFA -ot Int16  -te -168.5 15.0 -24.5 85.0 -ts 17280 8400 -r bilinear -srcnodata -9999 -dstnodata -9999 \

Flags: 

-of HFA   : Specify .img output file format. the most generic file format that supports image sizes > 2 GB (needed for our large mosaic sizes)

-ot Int16 : Output word size signed 16 bit - supports elevation values -10000 - 10000 meters

-te '-168.5 15.0 -24.5 85.0': Specify the bounding coordinates of the grid into which all of the input image 'tiles' will be resampled.

-ts 17280 8400: number of columns (X) and rows (Y) in the output grid: In combination with the grid definition (-te), determines the size of image pixels (in this case, 30 arcsecond @
3600 arc seconds in one degree): 

Columns (X):
   -168.5 - (-24.5)     = 144 degrees 
    144 degrees / 17280 = 0.0083333 degree pixel size (30 arcseconds or ~ 90 meters)

Rows (Y): 
    85 - 15             = 70 degrees
    70 / 8400           = 0.0833333 degree pixel size

-r bilinear: select bilinear resampling

-srcnodata -9999 -dstnodata -9999 : specify values for 'no data' in the input and output images

To execute the workflow: run scripts 1) through 3) sequentially.

File format note: The Tif file format only supports maximum image size of 2 GB (some platforms support 4 GB). Our image mosaics will exceed this limit, so we will use .IMG
format for the mosaiced images.

5/13 Note: Due to system load, I have not yet re-tested this workflow since changing
to the .img format. I will, ASAP. But this workflow has been successfully tested using .tif files, so it is a good example of the workflows to be used.
5/13 Note: I will most likely be replacing scripts 1) and 2) with virtual image tables, created using the GDAL 'gdalbuildvrt' utility. Using VRTs to represent the sub-mosaics
is more efficient as it eliminates the need for physical sub-mosaic files.

I am currently testing the VRT approach to make sure that it will work with the large mosaic file sizes encountered here.

When creating the ASTER GDEM image mosaics for the northeastern hemisphere, I noticed that the (sub) mosaics extended beyond the longitude boundaries implied
by the set of ASTER image tiles selected for the mosaic. This is easy to spot;
the mosaic included a large area of solid grey, with a single 1 degree ASTER
image tile embedded in the grey area (usually at the outer boundary).

Note that the dimensions and boundaries of these mosaics are defined by the
aggregate dimensions of the ASTER GDEM tiles included in the .vrt file for
the mosaic.

I checked the file list used to construct the .vrt for the 'oversized' mosaic;
the latitude/longitude coordinates embedded within each filename were consistent
with the intended mosaic boundaries.

Then, I checked the actual boundaries of each GDEM tile using gdalinfo. I discovered that five GDEM files actually covered areas (based on internal metadata) DIFFERENT than the areas implied by the file names:
Examples of the mismatch:

ASTGTM_N60E076_dem.tif ACTUALLY covers area centered on 91.5 long, 64.5 lat
ASTGTM_N61E075_dem.tif ACTUALLY covers area centered on 97.5 long, 65.5 lat
ASTGTM_N80E080_dem.tif ACTUALLY covers area centered on 91.5 long, 81.5 lat

This situation appears to occur in ASTER images for other parts of the Eastern Hemisphere; this problem delayed me in creating the baseline and derived data layer products that I had promised by June 4.

I need to determine which of the two ASTER data portals (NASA or JAPAN) supplied the erroneous ASTER files. If I can isolate the problem to one portal, then it will be easy to identify and replace the defective images with 'clean' images from the error-free repository. However, if BOTH data portals contain incorrectly tagged data tiles, then we will have to add an additional two steps to the mosaic production process: 1) obtain the image metadata for each tile with gdalinfo utility. 2) review the metadata for all ASTER images, 3) identify 'mistagged' ASTER images 4) replace the images, if possible with correct image tiles.

It is possible that no replacements are available for the mistagged files, and that the mosaics will contain some gaps.

Next step: Attempt to fix these problems within my two test mosaic areas.

Reviewing the 'NorthWestEast' ASTER GDEM component of the terrain mosaic,
I encountered one more mislabeled tile:

ASTGTM_N63E061_dem.tif

I replaced it with a 'correct' version from ASTER Japan portal, and recreated the image mosaic.

When I inspected the ASTER GDEM component of the North Eastern Hemisphere mosaic component,
I detected a series of defects in the image corresponding to a series of ASTER image tiles.
These appear as square or rectangular discontinuities extending throughout the lower half
of the image.

Reviewing the metadata for the tiles corresponding to the 'dropout' locations, I note that
many of them have location metadata that does not match the locations in the filenames.
Others have correct metadata, but do not match their surroundings.

To resolve it, I downloaded ~130 replacement files from the Aster Japan portal. These are
stored on /jupiter in the folder /data/organisms\rcr\ReplaceDefectiveTilesJuly2011/AsterDataJuly4.

I copied the replacement files to /home/reeves/EandO/asterGdem AFTER I copied the files being replaced to the (vulcan) subfolder ./IncorrectTilesJuly2011/repalcedJuly4.

Am re-running the mosaic building script, will review the modified image when it is done.

A global 90m product is now available, though the individual pieces have not been mosaiced together due to file size. Metadata that is correct for end-use distribution has not yet been prepared.

Persisting issues:

• Data were not pre-processed other than to shift the grids (SRTM and ASTER GDEM2) to align properly. No smoothing has been performed, and we are working on figuring out how to do this.
• Striping is apparent in the ASTER GDEM2 portion of the dataset (above N60), especially over Russia. We are discussing possible null pixel-filling procedures
• There is definite noise in the dataset, with pixel values as low as -200 (and probably lower). Perhaps smoothing the dataset will help, though more steps might need to be taken to
  fix the issue.