Reprojecting DTM from Graticule to Grid

We recently needed to process a digital terrain model (DTM) that was provided in an unprojected grid – ie in longitude and latitude and use in TrueViewVisuals. This post explains our process.

Data Validation

We use GDAL tools (part of OSGeo) to process raw data and also QGIS for desktop GUI with GRASS tools. For example, with new data the first thing we do is interrogate the data to validate the format and coordinate system using “gdalinfo” from the OSGeo shell.

gdalinfo w001001.adf

This returns metedata for the file and confirms we are working with WGS 84 data, for example:

Driver: AIG/Arc/Info Binary Grid
Files: .
Size is 5622, 2795
Coordinate System is:
GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (144.897477897848006,-38.324069951452998)
Pixel Size = (0.000045043421000,-0.000045043421000)
Corner Coordinates:
Upper Left ( 144.8974779, -38.3240700) (144d53'50.92"E, 38d19'26.65"S)
Lower Left ( 144.8974779, -38.4499663) (144d53'50.92"E, 38d26'59.88"S)
Upper Right ( 145.1507120, -38.3240700) (145d 9' 2.56"E, 38d19'26.65"S)
Lower Right ( 145.1507120, -38.4499663) (145d 9' 2.56"E, 38d26'59.88"S)
Center ( 145.0240950, -38.3870181) (145d 1'26.74"E, 38d23'13.27"S)
Band 1 Block=256x16 Type=Float32, ColorInterp=Undefined
Min=-0.980 Max=306.348
NoData Value=-3.4028234663852886e+38

To use the data in a regular grid requires not just a transformation (reprojection) but also a realignment. To transfom the “gdalwarp” command is available, but using QGIS makes the projection parameters much easier to manage (and we will also use QGIS for the grid realignment).

QGIS – New Project

Create a new project in QGIS and use the Project/Properties menu to apply a coordinate system. For this project I am using EPSG:28355 (GDA94 / MGA Zone 55) for data near Melbourne in Australia.

Now use the Layer/Add Layer/Add Raster Layer… menu to select the raster DTM file. You should see the raster stylised and appear as a trapezoid with non-vertical sides because of the projection (if you see the raster as a rectangle, check the projection in the project properties as above).

QGIS (GDAL) – Transform

We now have a few steps to turn the data into a regular grid (ie with 100m grid cell values).

  • Use the Raster/Projections/Warp (Reproject) menu and set the source CS (coordinate system) to ESPG:4326 WGS 84 and the target CS to the required projection.
    A resulting image sample is shown below with white values for the no-data (this example is a peninsula with sea on either side).

QGIS – Save and Export

  • From the ‘reprojected’ layer, use the right-click Export/Save As… menu to bring up the output options.
  • Within this dialog box of settings with region overrides to create a regular grid. Set regular extent boundaries and a grid resolution as per the example image below.
  • Finally, use the Raster/Conversion/Traslate (Convert Format)… menu to Export to an .ASC file. Looking at the header of this file, we see the regular grid that we expect.

ncols 1120
nrows 720
xllcorner 316200.000000000000
yllcorner 5742200.000000000000
cellsize 20.000000000000