Differences

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

--- nsw_dems [2021/02/10 14:03] – bushwalking
+++ nsw_dems [2023/10/08 21:47] – bushwalking
@@ Line 1: / Line 1: @@
 ====== Managing large DEMs ======
-===== Downloading data =====
+===== 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.
@@ Line 22: / Line 26: @@
 |8828_goulburn|8928_moss_vale|9028_kiama| | | | |
 |8827_braidwood|8927_ulladulla|9027_jervis_bay| | | | |
-===== Pre-processing data =====
-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 lowest -input_file_list "$targetFolder\input.txt" "$targetFolder\temp.vrt"
-</file>
-It is possible to then build a larger virtual raster from the individual 1:100k virtual rasters.
 ===== Loading data =====
@@ Line 169: / Line 38: @@
 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 })
 </code>
+===== 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 <code>-tap</code> flag set
+    * and possibly the <code>-ps</code> flag as well
+You can then use this merged raster for operations around the zone boundaries.