DHM25_in_SAGA_GIS

The digital height model DHM25 is a set of data representing the 3D form of the earth’s surface without vegetation and buildings. It is essentially based on the Swiss National Map 1:25 000.

Years ago, when Swisstopo was called Bundesamt für Landestopographie, they published this data in their home-made format called MMBLT, stored in files with an .mlt extension, who look like this from the inside:

MMBLT

A classical geodata-archiving problem arises: no software is able to read this format directly any more. There used to be r.in.swisstopo : a plugin for Grass GIS able to do the job, but that was Grass 6. Installation in Grass 7 simply doesn’t seem to work.

In 2005, Daniel Roethlisberger published a Perl script capable of converting MMBLT to PNG images and POV-Ray. Unfortunately, the script does not add world files (.pgw) to the PNG images. When I’ve tried to import them to qGIS or SAGA, all png tiles are pasted to the same x=0, y=0 origin, i.e. over each other.

So I’ve adapted the script to produce another format, standard in 2015:  an ARC/ESRI GRID in ASCII format. This data format looks like this:

ASC

Basically, as you can see, only the preamble changes. The file gets an .asc extension. It can then be simply dragged and dropped to SAGA GIS or qGIS for further treatment. In qGIS, when importing several .asc files, I particularly recommend using the tool Raster > Miscellaneous > Build virtual raster. The resulting .vrt file can also be dragged and dropped into SAGA for further treatment, like 3D visualisation or conversion of the elevation model into contour lines.

Here is the Perl code you need for MMBLT to ARC/ESRI GRID conversion

#!/usr/bin/env perl
# MMBLT to ASCII GRID - Convert SwissTopo/L+T MMBL(T) DHM25 elevation data
# http://ourednik.info
# by Andre Ourednik <andre.ourednik@unine.ch> Creative Commons BY 2015
# based on :
# MMBLT to POV-Ray by Daniel Roethlisberger <daniel@roe.ch> Copyright (C) 2004-2005
# http://www.roe.ch/DHM25
# 
#  
# Redistribution and use, with or without modification, are permitted
# provided that the following conditions are met:
# 1. Redistributions must retain the above copyright notice, this list of
#    conditions and the following disclaimer.
# 2. The name of the author may not be used to endorse or promote products
#    derived from this software without specific prior written permission.
# 
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
# IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
# OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
# NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
# THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
# 
# $Id: mmblt2ESRI.pl,v 1.3 2015/11/28 11:25:53 ourednik Exp $
# $Id: mmblt2pov.pl,v 1.3 2005/05/10 11:25:53 roe Exp $
use Getopt::Std;
# use GD;
use strict;

use vars qw{ $opt_m $opt_p $opt_t $opt_o $opt_q $opt_h };
getopts("m:p:t:o:qh");
sub usage
{
    die "$0 [OPTIONS...] [INFILE]

Read a DHM25 matrix model in MMBLT/MMBL format from INFILE or stdin, and
convert it to a more universally useful PNG height map and POV-Ray object.

	-o file   Write ASCII GRID to file
	-q        Be quiet
	-h        This cruft

If -o is not given, no output will be written, but the
MMBLT header and data will still be parsed and checked.

To transform mm1091.mlt into corresponding 1091_hm.png and 1091.inc, use:
	$0 -o 1091 mm1091.mlt
";
}
usage() if($opt_h);

my $heightmap_file = ($opt_o ? $opt_o .'.asc' : '');

my $in_file = shift() || '-';

sub coord
{
	my $c = shift() || die;
	$c =~ s/\.0+$//;
	$c =~ s/([0-9]{2,3})([0-9]{3})/$1\/$2/;
	return $c;
}

if($in_file eq '-')
{
	*INFILE = \*STDIN;
}
else
{
	open(INFILE, "<$in_file") or die "Cannot open $in_file for reading!";
}

#NEWHEADER
#--------------------------------------------------------------------------------
#DHM25-MATRIXMODELL                             (c)BUNDESAMT F. LANDESTOPOGRAPHIE
#--------------------------------------------------------------------------------
#NORD-WEST ECKE     [M]  XXXX00.0  XXXX00.0     ERSTER HOEHENWERT
#SUED-OST ECKE      [M]  XXXX00.0  XXXX00.0     LETZTER HOEHENWERT
#MASCHENWEITE WE/NS [M]      25.0      25.0
#MATRIXDIMENSIONEN WE/NS    XXX       XXX       TOTAL    XXXXXX MATRIXPUNKTE
#HOEHENBEREICH     [DM]    XXXX      XXXX       (6 CHARACTER PRO HOEHENWERT)
#--------------------------------------------------------------------------------
#FORMAT                   ASCII                 L+T-FORMAT DHM25-MATRIXMODELL
#RECORDLAENGE(CHAR.)       XXXX                 XXX HOEHENWERTE PRO RECORD
#--------------------------------------------------------------------------------
#ENDHEADER
#  XXXX  XXXX  XXXX  [...]

my $header;
while(($_ = <INFILE>) && ($_ !~ /^NEWHEADER$/)) { }
die "No MMBLT header found!\n" unless($_);

$header = $_; # <- NEWHEADER
while(($_ = <INFILE>) && ($_ !~ /^ENDHEADER$/)) { $header .= $_; }
$header .= $_; # <- ENDHEADER

# FORMAT                   ASCII                 L+T-FORMAT DHM25-MATRIXMODELL
$header =~ /FORMAT +ASCII +L\+T-FORMAT DHM25-MATRIXMODELL/ or die "Not a DHM25 ASCII MMBLT file!";

# NORD-WEST ECKE     [M]  XXXX00.0  XXXX00.0     ERSTER HOEHENWERT
$header =~ /NORD-WEST ECKE +\[M\] +([0-9.]+) +([0-9.]+)/ or die "No N/W edge in MMBLT header!";
my $edge_x1 = $1;
my $edge_y1 = $2;

# SUED-OST ECKE      [M]  XXXX00.0  XXXX00.0     LETZTER HOEHENWERT
$header =~ /SUED-OST ECKE +\[M\] +([0-9.]+) +([0-9.]+)/ or die "No S/E edge in MMBLT header!";
my $edge_x2 = $1;
my $edge_y2 = $2;

# MASCHENWEITE WE/NS [M]      25.0      25.0
$header =~ /MASCHENWEITE WE\/NS \[M\] +([0-9.]+) +([0-9.]+)/ or die "No grid size in MMBLT header!";
my $grid_x = $1;
my $grid_y = $2;

# MATRIXDIMENSIONEN WE/NS    XXX       XXX       TOTAL    XXXXXX MATRIXPUNKTE
$header =~ /MATRIXDIMENSIONEN WE\/NS +([0-9]+) +([0-9]+) +TOTAL +([0-9]+) MATRIXPUNKTE/ or die "No dimensions in MMBLT header!";
my $dim_x = $1;
my $dim_y = $2;
my $total = $3;

# HOEHENBEREICH     [DM]    XXXX      XXXX       (6 CHARACTER PRO HOEHENWERT)
$header =~ /HOEHENBEREICH +\[DM\] +([0-9]+) +([0-9]+)/ or die "No height range in MMBLT header!";
my $h_min = $1;
my $h_max = $2;

die "Inconsistent header data ($dim_x * $dim_y != $total)!" if($dim_x * $dim_y != $total);

my $coord_x1 = coord($edge_x1);
my $coord_y1 = coord($edge_y1);
my $coord_x2 = coord($edge_x2);
my $coord_y2 = coord($edge_y2);

my $factor_x = ($dim_x - 1) * $grid_x / 1000;
my $factor_y = ($dim_y - 1) * $grid_y / 1000;
my $factor_h = ($h_max - $h_min) / 10000;
my $offset_x = $edge_x1 / 1000;
my $offset_y = $edge_y2 / 1000;
my $offset_h = $h_min / 10000;

print <<EOF unless($opt_q);
MMBLT Source File   : $in_file
Data Points in Model: $dim_x x $dim_y = $total
Data Grid Size   (m): $grid_x x $grid_y
Height Range    (dm): $h_min .. $h_max
Stretch Factors (km): $factor_x / $factor_y / $factor_h
N/W Edge Coordinates: $coord_x1//$coord_y1
S/E Edge Coordinates: $coord_x2//$coord_y2
EOF


my $asc = "ncols\t" . $dim_x . "\n" 
	. "nrows\t" . $dim_y . "\n" 
	. "xllcorner\t" . $edge_x1 . "\n" 
	. "yllcorner\t" . $edge_y2 . "\n" 
	. "cellsize\t" . $grid_x ;

my $line;
for(my $y = 0; $y < $dim_y; $y++)
{
	$asc = $asc . "\n";
	for(my $x = 0; $x < $dim_x; $x++)
	{
		$line =~ s/^[\s\n\r]+//;
		$line = <INFILE> unless($line);
		die "Premature end of data!" unless($line);
		$line =~ /([0-9]+)/ or die "Illegal data!";
		$line = $';
		my $h = $1;
		die "Height value $h out of range $h_min..$h_max!" if($h < $h_min || $h > $h_max);
		$asc = $asc . $h;
		$asc = $asc . "\t";
	}
}

if($in_file ne '-')
{
	close(INFILE);
}

if($heightmap_file)
{
	open(OUTFILE, ">$heightmap_file") or die "Cannot open $heightmap_file for writing!";
	print OUTFILE $asc;
	close(OUTFILE);
}Code language: PHP (php)

To execute this code, copy it to your favourite text editor, save it as a file (named for example “mmblt2ESRI.pl-1.3”) to a directory of your choice, make it executable, and execute it like this:

perl mmblt2ESRI.pl-1.3 -o your_asc_output_file your_mlt_input_file.mlt

In order to convert all *.mlt files in a directory use this code (that will only work if all your .mlt files as well as your mmblt2ESRI.pl-1.3 script are in the same directory, and if this directory is also your current working directory) :

#!/bin/bash
for filename in *.mlt; do
perl mmblt2ESRI.pl-1.3 -o $(echo $filename | cut -f 1 -d '.') $filename
done

Of course, this code, too, can be saved in a batch file (called for example convert_to_asc.sh) and called from the commande line:

sudo chmod +x convert_to_asc.sh
./convert_to_asc.sh 

P.S.

Just for fun, this is what I did with the elevation model, using SAGA, qGIS and the qGIS2Leaf extension: A model of the Swiss territory with sea level at 1000m, as described in the Gérimont novels of Stéphane Bovon: Here a scrrenshot of the Napf Island:

NapfInsel1t
Napfinsel2t

Leave a comment

Your email address will not be published. Required fields are marked *