Managing large DEMs

(Note that this page can't be edited due to mod_security issues)

Intro

While data can be downloaded in an ad hoc manner, if you are regularly processing DEMs, it is better to have the DEM tiles already downloaded.

Once downloaded, pre-process the tiles into GeoTIFFs, build a virtual raster for each 1:100k grid square. You can then build a virtual raster of virtual rasters!

Downloading data

Download tiles by 1:100k map area, which is 0.5 x 0.5 degree squares. Each 1:100k map area ranges from around 2GB to 6GB of data, depending on the number of 2m DEMs vs 1m DEMs, and other factors.

The Sydney basin and Blue Mountains is around 50GB all up.

For example, to request data from ELVIS for the Katoomba 1:100k map area (bounded by -33.5,150,-34,150.5), see the screenshot below, and Select All DEMs:

For reference, below are the 1:100k map areas around Sydney:

8833_gulgong	8933_merriwa	9033_muswellbrook	9133_camberwell	9233_dungog	9333_buladelah	9433_forster
8832_mudgee	8932_mt_pomany	9032_howes_valley	9132_cessnock	9232_newcastle	9332_port_stephens
8831_bathurst	8931_wallerawang	9031_st_albans	9131_gosford	9231_lake_macquarie
8830_oberon	8930_katoomba	9030_penrith	9130_sydney
8829_taralga	8929_burragorang	9029_wollongong	9129_port_hacking
8828_goulburn	8928_moss_vale	9028_kiama
8827_braidwood	8927_ulladulla	9027_jervis_bay

And the boundaries:

Map	West	South	North	East
8833 Gulgong	149.5	-32.5	-32	150
8832 Mudgee	149.5	-33	-32.5	150
8831 Bathurst	149.5	-33.5	-33	150
8830 Oberon	149.5	-34	-33.5	150
8829 Taralga	149.5	-34.5	-34	150
8828 Goulburn	149.5	-35	-34.5	150
8827 Braidwood	149.5	-35.5	-35	150
8933 Merriwa	150	-32.5	-32	150.5
8932 Mt Pomany	150	-33	-32.5	150.5
8931 Wallerawang	150	-33.5	-33	150.5
8930 Katoomba	150	-34	-33.5	150.5
8929 Burragorang	150	-34.5	-34	150.5
8928 Moss Vale	150	-35	-34.5	150.5
8927 Ulladulla	150	-35.5	-35	150.5
9033 Muswellbrook	150.5	-32.5	-32	151
9032 Howes Valley	150.5	-33	-32.5	151
9031 St Albans	150.5	-33.5	-33	151
9030 Penrith	150.5	-34	-33.5	151
9029 Wollongong	150.5	-34.5	-34	151
9028 Kiama	150.5	-35	-34.5	151
9027 Jervis Bay	150.5	-35.5	-35	151
9133 Camberwell	151	-32.5	-32	151.5
9132 Cessnock	151	-33	-32.5	151.5
9131 Gosford	151	-33.5	-33	151.5
9130 Sydney	151	-34	-33.5	151.5
9129 Port Hacking	151	-34.5	-34	151.5
9233 Dungog	151.5	-32.5	-32	152
9232 Newcastle	151.5	-33	-32.5	152
9231 Lake Macquarie	151.5	-33.5	-33	152
9333 Buladelah	152	-32.5	-32	152.5
9332 Port Stephens	152	-33	-32.5	152.5
9433 Forster	152.5	-32.5	-32	153
8729 Crookwell	149	-34.5	-34	149.5
8728 Gunning	149	-35	-34.5	149.5

Pre-processing data

As of late 2023, data extracts are provided in a single zipped file with folders for each resolution (1m/2m/5m), with .TIF files in each folder. If you have data prior to late 2023 that you want to reprocess, see Pre-processing data (old) below.

Extract all files into a single folder.

These .TIF files have a file format like: Penrith201904-LID1-AHD_2806276_56_0002_0002_1m.tif

The '56' refers to the Zone. Different zones can't be mixed when creating a virtual raster, so any files with a different zone will need to be deleted or moved to another directory.

If the GDAL binaries are in your path, you can then run the following two commands in the folder at the (Windows) command line:

dir /b *.tif > input.txt
gdalbuildvrt.exe -resolution user -tr 2 2 -input_file_list input.txt area.vrt

MacOS and Linux users will need to convert these commands for their own use.

The area.vrt file is now a virtual raster of all of the other files.

It is possible to then build a larger virtual raster from the individual 1:100k virtual rasters. Note that larger rasters need to have the same projection. Also keep in mind that not all data is at 2m resolution.

Loading data

Loading up a large virtual raster into QGIS can be very slow, as can manipulating it. However, you can quickly load up a smaller section of the map using the following steps:

create the boundary of your desired area using the polygon tool at https://maps.ozultimate.com/
download the polygon using the Export drawn data to GeoJSON function
load up the Python console (Plugins→Python Console or Ctrl+Alt+P on Windows) in QGIS
run the following command (replace file locations with your own - the INPUT should be your virtual raster, and the MASK should be your GeoJSON boundary polygon)

resultClip = processing.runAndLoadResults("gdal:cliprasterbymasklayer", { 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : 'E:/geodata/8930_katoomba/area.vrt', 'KEEP_RESOLUTION' : False, 'MASK' : 'C:/Users/Downloads/data.geojson', 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : None, 'TARGET_CRS' : None, 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None })

If you don't want to use Python, you can use the Clip Raster by Mask Layer tool from the Processing Toolbox.

Multiple projections

If you are trying to process data near zone boundaries, you may run into issues where the tiles are of multiple projections. Virtual rasters run into numerous issues in this scenario. Even if you translate the data from one projection to the other, you will hit issues with pixel sizes, and pixel starting points.

The best approach so far seems to be to

clip each raster to a modest size
reproject one file to the other CRS (projection)
merge the two files, with
- the un-reprojected one first
- the -tap flag set
- and possibly the -ps flag as well

You can then use this merged raster for operations around the zone boundaries.

Pre-processing data (old)

This script was designed when the DEM files were delivered as doubly zipped .ASC files. As of late 2023, they come as .TIF files, albeit slightly larger files than the process below produced. This means that the steps below are now largely redundant. They are left in place in case people have old datasets that they want to process/reprocess. The current (2024) process is outlined above.

The following may be useful for Windows users.

Below is a Windows Powershell script that will

move any old DEMs and DEMs from a different zone (you can't mix zones in a virtual raster) to an archive sub-folder
extract the raw files from the remaining zip files
convert all of the .ASC files to GeoTIFF
move the old zip files to a current sub-folder
zip the current and archive sub-folders to temp.zip
create a virtual raster (.vrt) file of all of the GeoTiffs

You will need to replace the Environment variables (Path, GDAL_PATH, GDAL_DRIVER_PATH, PROJ_LIB) in the script with your own - see the lines starting with $Env.

Usage is:

.\buildvrt.ps1 <zipFile> <targetFolder>

eg

.\buildvrt.ps1 E:\geodata_raw\DATA_11279.zip E:\geodata\kiama

The .\ is required if you are running the script from your current directory.

buildvrt.ps1

$zipFile=$args[0]
$targetFolder=$args[1]
 
# Unzip files from all subdirectories to new folder
Expand-Archive -LiteralPath $zipFile -DestinationPath $targetFolder
Get-ChildItem -Path "$targetFolder\*.zip" -Recurse | Move-Item -Destination $targetFolder
Get-ChildItem -Path $targetFolder -Directory | Remove-Item -Recurse
 
# Create hash of zip files, by name (location, date)
$zipFileList = @{}
Get-ChildItem -Path "$targetFolder\*.zip" -Name | ForEach-Object {
	$zipFileList.add($_, @{}) 
	$_ -match '(\d{7})_(\d{2})'
	$zipFileList[$_]['location'] = $Matches.1
	$zipFileList[$_]['zone'] = $Matches.2
	$_ -match '[^_\d](\d{6})'
	$zipFileList[$_]['date'] = $Matches.1
}
# $zipFileList | ConvertTo-Json
 
# Create hash of location (date, name)
$locationList = @{}
$zoneCount = @{}
$zipFileList.keys | ForEach-Object {
	if($locationList[$zipFileList[$_]['location']]) {
		$t = $locationList[$zipFileList[$_]['location']]
		$t.add($zipFileList[$_]['date'],$_)
	} else {
		$t = @{}
		$t.add($zipFileList[$_]['date'],$_)
		$locationList.add($zipFileList[$_]['location'],$t)
	}
	if($zoneCount[$zipFileList[$_]['zone']]) {
		$zoneCount[$zipFileList[$_]['zone']]++
	} else {
		$zoneCount[$zipFileList[$_]['zone']] = 1
	}
}
# $locationList | ConvertTo-Json
# $zoneCount | ConvertTo-Json
 
# Create archive folder
$archiveFolder = "$targetFolder\archive"
If(!(test-path $archiveFolder))
{
      New-Item -ItemType Directory -Force -Path $archiveFolder
}
 
# Sort each location by date desc, and move old files to /archive
$locationList.keys | ForEach-Object {
	$i=0
	$locationList[$_].GetEnumerator() | sort key -des | ForEach-Object {
		if ($i -eq 0) {
			$i++
			return}
		else {
			#$_ | ConvertTo-Json
			$s = $_.Value
			Move-Item -Path "$targetFolder\$s" -Destination "$targetFolder\archive\$s"
		}
	}
}
 
# You can't build a VRT with files from a different projection, so
# delete files from outside main zone
# This could probably be altered to include a step to reproject those files
$mainZone = ''
$zoneCount.GetEnumerator() | sort value -des | select -first 1 | ForEach-Object {
	$mainZone = $_.Name
}
 
$zipFileList.keys | ForEach-Object {
	if($zipFileList[$_]['zone'] -ne $mainZone) {
		Remove-Item -Path "$targetFolder\$_"
	}
}
 
Get-ChildItem -Path "$targetFolder\*.zip" | Expand-Archive -DestinationPath $targetFolder
$Env:Path += ";C:\OSGEO4~1\apps\Python27\Scripts;C:\OSGEO4~1\bin"
$Env:GDAL_DATA="C:\OSGEO4~1\share\gdal"
$Env:GDAL_DRIVER_PATH="C:\OSGEO4~1\bin\gdalplugins"
$Env:PROJ_LIB="C:\OSGEO4~1\share\proj"
 
Get-ChildItem -Path "$targetFolder\*.asc" | ForEach-Object {
	$srcFile = $_.FullName
	$destFile = $_.FullName -replace '.asc', '.tif'
	&"gdal_translate.exe" -of "GTiff" $srcFile $destFile
}
 
Get-ChildItem -Path "$targetFolder\*.asc" | Remove-Item -Recurse
Get-ChildItem -Path "$targetFolder\*.html" | Remove-Item -Recurse
Get-ChildItem -Path "$targetFolder\*.prj" | Remove-Item -Recurse
Get-ChildItem -Path "$targetFolder\*.xml" | Remove-Item -Recurse
 
# Create current folder
$currentFolder = "$targetFolder\current"
If(!(test-path $currentFolder))
{
      New-Item -ItemType Directory -Force -Path $currentFolder
}
 
# Move zip files to current folder
Move-Item -Path "$targetFolder\*.zip" -Destination "$currentFolder"
 
# Zip /archive & /current to new zip folder
Compress-Archive -Path "$targetFolder\current", "$targetFolder\archive" -DestinationPath "$targetFolder\temp.zip"
 
# Create 2m vrt
New-Item "$targetFolder\input.txt"
Get-ChildItem -Path "$targetFolder\*.tif" | Add-Content "$targetFolder\input.txt"
&"gdalbuildvrt.exe" -resolution user -tr 2 2 -input_file_list "$targetFolder\input.txt" "$targetFolder\temp.vrt"

You can rename temp.vrt to anything you like, but it needs to stay relative to the GeoTIFF files.