[GRASS-SVN] r41564 - in grass-addons/raster: . r.in.swisstopo

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Mar 26 13:19:50 EDT 2010


Author: Juergen
Date: 2010-03-26 13:19:49 -0400 (Fri, 26 Mar 2010)
New Revision: 41564

Added:
   grass-addons/raster/r.in.swisstopo/
   grass-addons/raster/r.in.swisstopo/Makefile
   grass-addons/raster/r.in.swisstopo/description.html
   grass-addons/raster/r.in.swisstopo/r.in.swisstopo
Log:
script for importing swisstopo data

Added: grass-addons/raster/r.in.swisstopo/Makefile
===================================================================
--- grass-addons/raster/r.in.swisstopo/Makefile	                        (rev 0)
+++ grass-addons/raster/r.in.swisstopo/Makefile	2010-03-26 17:19:49 UTC (rev 41564)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.in.swisstopo
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

Added: grass-addons/raster/r.in.swisstopo/description.html
===================================================================
--- grass-addons/raster/r.in.swisstopo/description.html	                        (rev 0)
+++ grass-addons/raster/r.in.swisstopo/description.html	2010-03-26 17:19:49 UTC (rev 41564)
@@ -0,0 +1,104 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS manual: r.in.swisstopo</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>r.in.swisstopo</b></em>  - Imports a digital elevation model from swisstopo into a Grass raster map and interpolates it into the current mapset's resolution. The name of the raster map is defined by the user input. A further raster layer (with the appendix '_origres') is automatically created and contains the input data in its original resolution. Resultant raster maps will contain elevation data in meters.
+<h2>KEYWORDS</h2>
+swisstopo, DHM, DEM, Raster, Import
+<h2>SYNOPSIS</h2>
+<b>r.in.swisstopo</b><br>
+<b>r.in.swisstopo help</b><br>
+<b>r.in.swisstopo</b> <b>input</b>=<em>input</em> <b>output</b>=<em>raster</em>  [<b>method</b>=<em>string</em>]   [<b>type</b>=<em>string</em>]   [<b>zrange</b>=<em>min,max</em>]   [<b>percent</b>=<em>integer</em>]   [<b>method_resamp</b>=<em>string</em>]   [--<b>overwrite</b>]  [--<b>verbose</b>]  [--<b>quiet</b>] 
+
+<h3>Flags:</h3>
+<DL>
+<DT><b>--overwrite</b></DT>
+<DD>Force overwrite of output files</DD>
+<DT><b>--verbose</b></DT>
+<DD>Verbose module output</DD>
+<DT><b>--quiet</b></DT>
+<DD>Quiet module output</DD>
+</DL>
+
+<h3>Parameters:</h3>
+<DL>
+<DT><b>input</b>=<em>input</em></DT>
+<DD>Swisstopo file containing DEM data (*.mlt, *.mbl or *.xyz)</DD>
+
+<DT><b>output</b>=<em>raster</em></DT>
+<DD>Name of the raster layer, that will be created out of the imported data</DD>
+
+<DT><b>method</b>=<em>string</em></DT>
+<DD>Statistic to use for raster values</DD>
+<DD>Options: <em>n,min,max,range,sum,mean,stddev,variance,coeff_var</em></DD>
+<DD>Default: <em>mean</em></DD>
+
+<DT><b>type</b>=<em>string</em></DT>
+<DD>Storage type for resultant raster map</DD>
+<DD>Options: <em>CELL,FCELL,DCELL</em></DD>
+<DD>Default: <em>FCELL</em></DD>
+
+<DT><b>zrange</b>=<em>min,max</em></DT>
+<DD>Filter range for z data (min,max)</DD>
+
+<DT><b>percent</b>=<em>integer</em></DT>
+<DD>Percentage of map to keep in memory</DD>
+<DD>Options: <em>1-100</em></DD>
+<DD>Default: <em>100</em></DD>
+
+<DT><b>method_resamp</b>=<em>string</em></DT>
+<DD>Interpolation method for interpolation of input data to current region's resolution (using r.resamp.interp)</DD>
+<DD>Options: <em>nearest,bilinear,bicubic</em></DD>
+<DD>Default: <em>bilinear</em></DD>
+
+<h2>Manual download and installation of the module (without svn)</h2>
+<p>You can find the current version of r.in.swisstopo <a href="www.">here</a><br>
+<br>
+Installation:<br>
+You might need to have root privileges for the following. If you do not have these, please ask your system administrator
+<br>
+<ul>
+<li>Simply copy the file 'r.in.swisstopo.sh' into your $GISBASE/scripts/ directory.<br>
+(To find out your $GISBASE directory, call the following command from within a GRASS terminal:<br>
+<i>env | grep $GISBASE</i>)<br>
+Or directly copy the file from within the GRASS terminal with something like:
+<i>cp r.in.swisstopo.sh $GISBASE/scripts/</i>
+<br></li>
+<li>Additionally, you might have to make the script executable with:<br>
+<i>chmod +x $GISBASE/scripts/r.in.swisstopo.sh</i>
+<br></li>
+<li>Finally, copy this file (description.html) into the $GISBASE/docs/html/ directory and name it r.in.swisstopo.sh
+<br></li>
+</ul>
+<br>
+That's it!
+</p>
+
+<h2>SEE ALSO</h2>
+<p>
+<a href="r.in.xyz.html">r.in.xyz</a><br>
+<a href="r.resamp.interp.html">r.resamp.interp</a><br>
+</p>
+
+<h2>Note</h2>
+<p>If you find any bugs or find the module helpful or have any questions, feel free to contact the author (Juergen.Hansmann["at"]erdw.ethz.ch)</p>
+
+<h2>Author</h2>
+<p>J&uuml;rgen Hansmann<br>
+<i>Department of Earth Sciences<br>
+Swiss Federal Institute of Technology (ETH)<br>
+Zurich<br>
+</i></p>
+<p><i>Last changed: $Date: 03/26/2010 $</i> </p> 
+
+</DL>
+</body>
+</html>

Added: grass-addons/raster/r.in.swisstopo/r.in.swisstopo
===================================================================
--- grass-addons/raster/r.in.swisstopo/r.in.swisstopo	                        (rev 0)
+++ grass-addons/raster/r.in.swisstopo/r.in.swisstopo	2010-03-26 17:19:49 UTC (rev 41564)
@@ -0,0 +1,276 @@
+#!/bin/bash
+
+############################################################################
+#
+# MODULE:       r.in.swisstopo
+# AUTHOR(S):    Juergen Hansmann, 2010
+# COPYRIGHT:    (C) 2010 GRASS Development Team/Juergen Hansmann
+#
+# imports digital elevation models from swisstopo into a GRASS raster map
+# using r.in.xyz and r.resamp.interp
+#
+#
+#  This program is free software; you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation; either version 2 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#############################################################################/
+#%Module
+#% description: Imports a digital elevation model from swisstopo into a Grass raster map and interpolates it into the current mapset's resolution. The name of the raster map is defined by the user input. A further raster layer (with the appendix '_origres') is automatically created and contains the input data in its original resolution. Resultant raster maps will contain elevation data in meters.
+#% keywords: swisstopo, DHM, DEM, Raster, Import
+#%End
+#%Option
+#% key: input
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: input
+#% description: Swisstopo file containing DEM data (*.mlt, *.mbl or *.xyz)
+#% gisprompt: old_file,file,input
+#%End
+#%Option
+#% key: output
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: raster
+#% description: Name of the raster layer, that will be created out of the imported data
+#% gisprompt: new,cell,raster
+#%End
+#%Option
+#% key: method
+#% type: string
+#% required: no
+#% multiple: no
+#% options: n,min,max,range,sum,mean,stddev,variance,coeff_var
+#% description: Statistic to use for raster values
+#% answer: mean
+#%End
+#%Option
+#% key: type
+#% type: string
+#% required: no
+#% multiple: no
+#% options: CELL,FCELL,DCELL
+#% description: Storage type for resultant raster map
+#% answer: FCELL
+#%End
+#%Option
+#% key: zrange
+#% type: double
+#% required: no
+#% multiple: no
+#% key_desc: min,max
+#% description: Filter range for z data (min,max)
+#%End
+#%Option
+#% key: percent
+#% type: integer
+#% required: no
+#% multiple: no
+#% options: 1-100
+#% description: Percentage of map to keep in memory
+#% answer: 100
+#%End
+#%Option
+#% key: method_resamp
+#% type: string
+#% required: no
+#% multiple: no
+#% options: nearest,bilinear,bicubic
+#% description: Interpolation method for interpolation of input data to current region's resolution (using r.resamp.interp)
+#% answer: bilinear
+#%End
+
+
+
+########## some basic checks ###################################
+
+# check if we are in the GRASS environment
+if [ -z "$GISBASE" ] ; then
+    echo "You must be in GRASS GIS to run this program." 1>&2
+    exit 1
+fi
+
+
+# check if we have awk
+if [ ! -x "`which awk`" ] ; then
+    g.message -e "awk required, please install awk or gawk first" 
+    exit 1
+fi
+
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    exec g.parser "$0" "$@"
+fi
+
+
+# check if input file exists
+if [ ! -e $GIS_OPT_INPUT ] ; then
+  g.message -e "Input file <$GIS_OPT_INPUT> not found!"
+  exit 1
+fi
+
+
+# check, if the input file has the correct type
+file_extension=`echo $GIS_OPT_INPUT | cut -d . -f 2`
+
+case "$file_extension" in
+	"mlt" | "MLT" | "mbl" | "MBL" | "xyz" | "XYZ" ) g.message "importing swisstopo data. Depending on the file size (and your computer) this may take a while..";;
+	*) g.message -e "Input file has wrong type!"
+	exit 1;;
+esac
+
+################################################################/
+
+
+# setup temporary ascii file for storage of converted data
+TMP="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
+  g.message -e "unable to create temporary files"
+  exit 1
+fi
+
+
+# setup temporary file for storage of region settings
+TMP2="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP2" ] ; then
+  g.message -e "unable to create temporary files"
+  exit 1
+fi
+
+
+# set environment so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+
+
+# Now do the hard work
+if [ "$file_extension" == "xyz" -o "$file_extension" == "XYZ"  ] ; then
+
+
+# if the swisstopo input file is in .xyz format already, the only thing, that needs to be done is to
+# remove the different number of blanks between the columns with only one blank each.
+  awk '{
+	print $1, $2, $3;
+	if (FNR==1){last_x=$1;last_y=$2};
+	x_res=$1-last_x;
+	last_x=$1;
+	if(x_res<0){x_res=-x_res};
+	if($2 != last_y){y_res=$2-last_y;last_y=$2};
+	if(y_res<0){y_res=-y_res};	
+       }
+  END{
+  print x_res, y_res >> "'$TMP2'";
+  }' $GIS_OPT_INPUT > $TMP
+
+  # get the input file's x-resolution (W/E resolution) and y-resolution (N/S resolution)
+  xres=`awk < $TMP2 '{print $1}'`
+  yres=`awk < $TMP2 '{print $2}'`
+
+
+else
+
+# if the swisstopo input file is in .mlt or .mbl format, a bit more work needs to be done
+  awk ' BEGIN{x_count=0; y_count=0;inc_count=0;start_read=0};
+    {
+	  if (FNR == 3)
+	    {split($1, a, "-");type=a[2]}	
+
+	  if ($1 == "NORD-WEST")
+	    {min_x=$4; max_y=$5}
+
+	  if ($1 == "SUED-OST")
+	    {max_x=$4; min_y=$5}
+
+	  if ($1 == "MASCHENWEITE")
+	    {x_res=$4; y_res=$5}
+
+	  if ($1 == "MATRIXDIMENSIONEN")
+	    {matrix_increment=$3;}
+
+	  if ($1 == "ENDHEADER")
+	    {start_read=FNR}
+
+	  if (FNR > start_read && start_read != 0)
+	    {
+
+	      for (i=1; i<=NF; i++)
+		  {
+		  inc_count=inc_count+1;
+		  x_utm=min_x+(x_count*x_res);
+		  y_utm=max_y-(y_count*y_res);
+		  print x_utm,y_utm,$i*0.1;
+		  x_count=x_count+1
+		  if (inc_count == matrix_increment)
+		    {inc_count=0; x_count=0; y_count=y_count+1};
+		  }
+	    }
+    }
+  END{
+  print min_x, max_x, min_y, max_y, x_res, y_res, type >> "'$TMP2'";
+  }' $GIS_OPT_INPUT > $TMP
+
+
+  # check, wether header information of input file declares the file to be a "MATRIXMODELL"
+  # if it is not: exit!
+  type=`awk < $TMP2 '{print $7}'`
+  if [ "$type" != "MATRIXMODELL" ] ; then
+    g.message -e "Something seems to be wrong with the input file format"
+    exit 1
+  fi
+
+
+  # get the input file's x-resolution (W/E resolution) and y-resolution (N/S resolution)
+  xres=`awk < $TMP2 '{print $5}'`
+  yres=`awk < $TMP2 '{print $6}'`
+
+
+fi
+
+
+################## preparing to import the data into the current mapset #########################
+
+# get current GRASS resolution
+eval `g.region -g`
+# now two variables exist: nsres and ewres, containing the resolution of the current mapset
+
+
+# set resolution to the input dataset resolution
+g.region nsres=$yres ewres=$xres
+
+# import input data with r.in.xyz
+g.message "invoking r.in.xyz"
+
+
+if [ -z "$GIS_OPT_ZRANGE" ] ; then
+  r.in.xyz input=$TMP output=${GIS_OPT_OUTPUT}_origres method=$GIS_OPT_METHOD type=$GIS_OPT_TYPE fs="space" x=1 y=2 z=3 percent=$GIS_OPT_PERCENT
+else
+  r.in.xyz input=$TMP output=${GIS_OPT_OUTPUT}_origres method=$GIS_OPT_METHOD type=$GIS_OPT_TYPE fs="space" x=1 y=2 z=3 zrange=$GIS_OPT_ZRANGE percent=$GIS_OPT_PERCENT
+fi
+
+
+# reset resolution to original resolution of the mapset
+g.region nsres=$nsres ewres=$ewres
+
+r.resamp.interp input=${GIS_OPT_OUTPUT}_origres out=$GIS_OPT_OUTPUT meth=$GIS_OPT_METHOD_RESAMP
+
+
+#################################################################################################/
+
+
+# cleaning up
+rm $TMP
+rm $TMP2
+
+
+g.message "Done!"
+


Property changes on: grass-addons/raster/r.in.swisstopo/r.in.swisstopo
___________________________________________________________________
Added: svn:executable
   + *



More information about the grass-commit mailing list