Differences

This shows you the differences between two versions of the page.

--- nsw_dems [2023/10/08 21:47] – bushwalking
+++ nsw_dems [2023/11/23 15:01] (current) – [Multiple projections] bushwalking
@@ Line 1: / Line 1: @@
 ====== Managing large DEMs ======
+(Note that this page can't be edited due to mod_security issues)
 ===== Intro =====
@@ Line 26: / Line 28: @@
 |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 (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. See below for new steps.
+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:
+<code>.\buildvrt.ps1 <zipFile> <targetFolder></code>
+eg <code>.\buildvrt.ps1 E:\geodata_raw\DATA_11279.zip E:\geodata\kiama</code>
+The .\ is required if you are running the script from your current directory.
+<file powershell 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"
+</file>
+You can rename temp.vrt to anything you like, but it needs to stay relative to the GeoTIFF 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.
+===== Pre-processing data (new) =====
+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.
+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:
+<code>
+dir /b *.tif > input.txt
+gdalbuildvrt.exe -resolution user -tr 2 2 -input_file_list input.txt area.vrt
+</code>
+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 =====
@@ Line 48: / Line 250: @@
   * merge the two files, with
     * the un-reprojected one first
-    * the <code>-tap</code> flag set
+    * the ''-tap'' flag set
-    * and possibly the <code>-ps</code> flag as well
+    * and possibly the ''-ps'' flag as well
 You can then use this merged raster for operations around the zone boundaries.