I am attempting to take a simple opentopo heightmap and carry out the following in QGIS 3.4.1

- Reproject to UTM 31N
- Set cell size resolution to 5mtr
- Clip to a shapefile polygon feature exactly 40960 x 40960mtr

So why when using following command in Python console:

import os
os.system(r'''gdalwarp -t_srs EPSG:32631 -tr 5.0 5.0 -of GTiff -cutline E:\Users\rosso\Arma\QGIS\UTM31N\Test23-11-18\shape.shp -crop_to_cutline "E:/Users/rosso/Arma/Terrain Heightmaps/Heightmaps Opentopo/output_srtm.asc" E:/Users/rosso/Arma/QGIS/UTM31N/24-11-18/converted2.tif''')

Do I get width and height not matching exactly to my shapefile feature I'm clipping to - (40960 x 40960) - which at 5mtr resolution size should be '8192' for both width and Height

enter image description here

If I break up the process into 3 stages:

Change CRM:

import os
os.system(r'''gdalwarp -t_srs EPSG:32631 -of GTiff "E:/Users/rosso/Arma/Terrain Heightmaps/Heightmaps Opentopo/output_srtm.asc" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/crm.tif''')

Change Resolution:

import os
os.system(r'''gdalwarp -tr 5.0 5.0 -of GTiff "E:/Users/rosso/Arma/QGIS/UTM31N/25-11-18/crm.tif" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/res.tif''')

Clip to Shapefile feature mask layer:

import os
os.system(r'''gdalwarp -of GTiff -cutline E:\Users\rosso\Arma\QGIS\UTM31N\Test23-11-18\shape.shp -crop_to_cutline "E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/res.tif" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/cut.tif''')

The width and height adhere to what I expect:

enter image description here

I would prefer the process to be completed in one command, rather than splitting them out as I want to reduce the amount of processing on the original raster.

  • So I progressed a little in combining the first 2 commands for CRM and Resolution - gdalwarp -t_srs EPSG:32631 -wo SOURCE_EXTRA=1000 -tr 5.0 5.0 -r cubic -of GTiff "E:/Users/rosso/Arma/Terrain Heightmaps/Heightmaps Opentopo/output_srtm.asc" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/crm.tif But I still can't seem to combine with the Clip to Shapefile command without it messing with width and height. – Rossoe Dec 1 '18 at 12:12

I discovered that due to an update to GDAL (2.4.0) using 'Clip Raster by Mask Layer' no longer works as I expected - see thread I posted here - https://github.com/OSGeo/gdal/issues/1176

So I switched to following command and fixed my issue:

import os
os.system(r'''gdalwarp -t_srs EPSG:32630 -wo SOURCE_EXTRA=1000 -tr 5.0 5.0 -srcnodata "-9999" -r cubic -of GTiff -te 597700.541507 5593677.822312 638660.541507 5634637.822312 "E:/Arma/Heightmaps/Opentopo/output_srtm.asc" E:/Arma/project/converted.tif''')

Also followed a few tips from this chap: https://medium.com/planet-stories/a-gentle-introduction-to-gdal-part-2-map-projections-gdalwarp-e05173bd710a

particularly re some additional warp parameters that improved raster quality!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.