[GRASS-SVN] r43455 - in grass-addons/raster: . r.out.kap_template

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Sep 13 20:14:03 EDT 2010


Author: hamish
Date: 2010-09-14 00:14:03 +0000 (Tue, 14 Sep 2010)
New Revision: 43455

Added:
   grass-addons/raster/r.out.kap_template/
   grass-addons/raster/r.out.kap_template/Makefile
   grass-addons/raster/r.out.kap_template/r.out.kap_template
Log:
+module to write a GeoTiff and metadata file in a way compatible with KAP/BSB encoders

Added: grass-addons/raster/r.out.kap_template/Makefile
===================================================================
--- grass-addons/raster/r.out.kap_template/Makefile	                        (rev 0)
+++ grass-addons/raster/r.out.kap_template/Makefile	2010-09-14 00:14:03 UTC (rev 43455)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.out.kap_template
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script


Property changes on: grass-addons/raster/r.out.kap_template/Makefile
___________________________________________________________________
Added: svn:mime-type
   + text/x-makefile
Added: svn:eol-style
   + native

Added: grass-addons/raster/r.out.kap_template/r.out.kap_template
===================================================================
--- grass-addons/raster/r.out.kap_template/r.out.kap_template	                        (rev 0)
+++ grass-addons/raster/r.out.kap_template/r.out.kap_template	2010-09-14 00:14:03 UTC (rev 43455)
@@ -0,0 +1,629 @@
+#!/bin/sh
+############################################################################
+#
+# MODULE:       r.out.kap_template
+#
+# AUTHOR(S):    H. Bowman
+#
+# PURPOSE:      Create a KAP file header and accompanying TIFF file
+#
+# COPYRIGHT:    (c) 2010 the OpenCPN Development Team
+#
+#               This program is free software under the GNU General Public
+#               License (>=v2). Read the file COPYING that comes with GRASS
+#               for details.
+#
+#		*** NOT FOR NAVIGATIONAL USE - USE AT YOUR OWN RISK ***
+#		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: Create a KAP file header and accompanying TIFF file.
+#% keywords: raster, export
+#%End
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of input map
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new_file,file,file
+#% label: Name for output file (extension automatically added)
+#% description: default: input map name with .hdr extention
+#% required : no
+#%end
+#%option
+#% key: bounding_poly
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Name of vector map containing single bounding box polygon
+#% required : no
+#%end
+#%option
+#% key: scale
+#% type: integer
+#% description: Map scale (if not given it will try to calculate one)
+#% required : no
+#%end
+#%option
+#% key: dpi
+#% type: integer
+#% description: Map pixels per inch (if not given it will try to calculate one)
+#% required : no
+#%end
+#%option
+#% key: group
+#% type: string
+#% gisprompt: old,group,group
+#% description: Name of georectification group (XY locations only)
+#% required : no
+#% guisection: Simple XY
+#%end
+#%option
+#% key: projection
+#% type: string
+#% description: Map projection (XY locations only)
+#% options: mercator,transverse mercator,lambert conformal conic,polyconic
+#% answer: mercator
+#% required : no
+#% guisection: Simple XY
+#%end
+#%option
+#% key: lat_0
+#% type: double
+#% description: Latitude for map projection (if applicable; XY locations only)
+#% required : no
+#% guisection: Simple XY
+#%end
+#%option
+#% key: lon_0
+#% type: double
+#% description: Longitude for map projection (if applicable; XY locations only)
+#% required : no
+#% guisection: Simple XY
+#%end
+#%option
+#% key: datum
+#% type: string
+#% description: Map datum (XY locations only)
+#% answer: WGS84
+#% required : no
+#% guisection: Simple XY
+#%end
+#%flag
+#% key: e
+#% description: Export the header only
+#%end
+#%flag
+#% key: z
+#% description: Crop away as much collar as possible from the edges of the image
+#%end
+#%flag
+#% key: b
+#% description: Bypass the potentially slow color reduction step
+#%end
+#%flag
+#% key: d
+#% description: Chart crosses the international dateline
+#%end
+
+
+if  [ -z "$GISBASE" ] ; then
+   echo "You must be in GRASS GIS to run this program." >&2
+   exit 1
+fi   
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+   exec g.parser "$0" "$@"
+fi
+
+
+g.message -i "This is EXPERIMENTAL software. NOT FOR NAVIGATIONAL USE."
+g.message message="
+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. Use at your own risk."
+
+
+MAP_NAME="$GIS_OPT_INPUT"
+
+if [ -n "$GIS_OPT_OUTPUT" ] ; then
+   OUTFILE="$GIS_OPT_OUTPUT"
+else
+   OUTFILE="$MAP_NAME"
+fi
+OUTFILE=`echo "$OUTFILE" | sed -e 's/\.hdr$//'`
+
+
+### check if outfile already exists and $GRASS_OVERWRITE ###
+### and check for .tif if -e flag was not given
+
+MAP_DESCRIPTION=`r.info -m "$MAP_NAME" | cut -f2- -d=`
+NU=`echo "$MAP_NAME" | tr '[:lower:]' '[:upper:]' | tr -d '[:punct:]'`
+  #??? #MAP_NAME=  CHART NUMBER?
+
+
+if [ -z "$GIS_OPT_DPI" ] && [ -z "$GIS_OPT_SCALE" ] ; then
+   g.message -e "Either dpi or scale must be given (both is ok too)."
+   exit 1
+fi
+
+
+PROJ_TYPE=`g.region -p | grep '^projection' | cut -f2 -d' '`
+# 0 : XY
+# 99 : other
+if [ "$PROJ_TYPE" -eq 0 ] ; then
+   # XY location
+   GROUP=`echo "$GIS_OPT_GROUP" | cut -f1 -d'@'`
+   # FIXME: support @othermapset
+   if [ -z "$GROUP" ] ; then
+      g.message -e "Georectifation group required for XY data source."
+      exit 1
+   fi
+   if [ -z "$GIS_OPT_PROJECTION" ] ; then
+      g.message -e "Map projection required for XY data source."
+      exit 1
+   fi
+   g.findfile element=group/$GROUP file=POINTS > /dev/null
+   if [ $? -ne 0 ] ; then
+      g.message -e "Could not find ground control POINTS file for this group."
+      exit 1
+   fi
+   if [ -z "$GIS_OPT_DPI" ] || [ -z "$GIS_OPT_SCALE" ] ; then
+      g.message -e "Both dpi and scale must be given for non-geographic use"
+      exit 1
+   fi
+   if [ "$GIS_FLAG_Z" -eq 1 ] ; then
+      g.message -e "XY maps must be exported at their full dimension"
+      exit 1
+   fi
+
+fi
+
+
+if [ "$GIS_FLAG_Z" -eq 1 ] ; then
+   # (was decollar_noaa_charts.sh)
+   #    http://lists.osgeo.org/pipermail/grass-user/2009-April/049723.html
+   ###simple standby alternative:
+   ###g.region zoom="$MAP_NAME"
+
+   # check that we are dealing with a paletted map
+   eval `r.info -t "$MAP_NAME"`
+   if [ "$datatype" != "CELL" ] ; then
+      echo "ERROR: only categorical maps can be processed for collar removal"
+   fi
+
+   # find cats in map
+   CATS=`r.category $MAP_NAME`
+
+   # find black and white category numbers
+   #   (salmon used in some NZ charts, repurpose as needed)
+   unset BLACK WHITE SALMON
+   for CAT in $CATS ; do
+     RULE=`r.what.color in="$MAP_NAME" value="$CAT"`
+     g.message -v "$RULE"
+     COLOR=`echo "$RULE" | cut -f2 -d' '`
+     if [ "$COLOR" = "0:0:0" ] ; then
+        BLACK="$CAT"
+     elif [ "$COLOR" = "255:255:255" ] ; then
+        WHITE="$CAT"
+     elif [ "$COLOR" = "219:73:150" ] ; then
+        SALMON="$CAT"
+     fi
+   done
+   g.message -v "black is $BLACK, white is $WHITE, salmon is $SALMON"
+
+   REMAINING=`echo $CATS | tr ' ' '\n' | grep -vw "$BLACK\|$WHITE\|$SALMON"`
+
+   #### setup temporary file
+   TMPFILE="`g.tempfile pid=$$`"
+   if [ $? -ne 0 ] || [ -z "$TMPFILE" ] ; then
+       g.message -e "Unable to create temporary file"
+       exit 1
+   fi
+
+   echo "$BLACK = NULL" > "$TMPFILE"
+   echo "$WHITE = NULL" >> "$TMPFILE"
+   echo "$SALMON = NULL" >> "$TMPFILE"
+   for CAT in $REMAINING ; do
+      echo "$CAT = $CAT" >> "$TMPFILE"
+   done
+
+
+   # remove stray single dots
+   r.neighbors in="$MAP_NAME" out="$MAP_NAME.despeckle_$$" method=mode
+
+   # make a reclass map, setting black and white to NULL
+   r.reclass in="$MAP_NAME.despeckle_$$" out="$MAP_NAME.filtered_$$" rules="$TMPFILE"
+
+   g.region zoom="$MAP_NAME.filtered_$$"
+   r.mapcalc "$MAP_NAME.cropped_$$ = $MAP_NAME"
+
+   #r.mapcalc "$MAP_NAME.cropped = $MAP_NAME.filtered"
+   
+   ## debug method for finding residuals
+      #r.mapcalc "$MAP_NAME.cover = 1"
+      #g.region rast="$MAP_NAME"
+      #r.mapcalc "$MAP_NAME.residual = if($MAP_NAME.cover == 1, null(), $MAP_NAME)"
+      #r.stats -c "$MAP_NAME.residual"
+      ## single 6,8 dot
+      #d.zoom to collar which should be empty
+      #r.univar $MAP_NAME
+      #r.out.xyz $MAP_NAME
+      #echo "symbol basic/star 20 -72.35262972 38.8881031 black blue" | d.graph -m
+      #echo "symbol basic/star 20 -74.51749113 39.27423606 black aqua" | d.graph -m
+      
+   #r.out.gdal in="$MAP_NAME" out="${MAP_NAME}.decollared.tif" type=Byte #...
+
+   # cleanup
+   g.remove rast="$MAP_NAME.filtered_$$" --quiet
+   g.remove rast="$MAP_NAME.despeckle_$$" --quiet
+   \rm -f "$TMPFILE"
+
+   # reset target input image to cropped version
+   MAP_NAME="$MAP_NAME.cropped_$$"
+fi
+
+eval `g.region -g`
+
+
+#### set up dpi and map scale terms
+
+if [ -n "$GIS_OPT_DPI" ] ; then
+   DPI="$GIS_OPT_DPI"
+else
+   # calc from map extent and scale
+   #dpi = rows * scale * i2m / ns_extent
+   eval `g.region -eg`
+   MAP_SCALE=`echo "$GIS_OPT_SCALE" | cut -f2 -d:`
+
+   DPI_NS=`echo "$rows $MAP_SCALE $ns_extent" | awk '{print $1 * $2 * 0.0254 / $3}'`
+   DPI_EW=`echo "$cols $MAP_SCALE $ew_extent" | awk '{print $1 * $2 * 0.0254 / $3}'`
+   # the two above should be the same.
+   DPI=`echo "$DPI_NS $DPI_EW" | awk '{print int(0.5 + ($1 + $2) / 2)}'`
+fi
+
+if [ -n "$GIS_OPT_SCALE" ] ; then
+   # fix in case the user gave 1:xxxxx
+   MAP_SCALE=`echo "$GIS_OPT_SCALE" | cut -f2 -d:`
+else
+   # calc from map extent and dpi
+   #SCALE_NS=$ns_extent/(($rows/$DPI)*$INCH2METER)
+   #SCALE_EW=$ew_extent/(($cols/$DPI)*$INCH2METER)
+   eval `g.region -eg`
+   DPI="$GIS_OPT_DPI"
+
+   SCALE_NS=`echo "$ns_extent $rows $DPI" | awk '{print $1/(($2/$3)*0.0254)}'`
+   SCALE_EW=`echo "$ew_extent $cols $DPI" | awk '{print $1/(($2/$3)*0.0254)}'`
+   # the two above should be the same.
+   MAP_SCALE=`echo "$SCALE_NS $SCALE_EW" | awk '{print int(0.5 + ($1 + $2) / 2)}'`
+fi
+
+
+#### set up projection terms
+if [ "$PROJ_TYPE" -eq 0 ] ; then
+   # XY location
+   PROJECTION=`echo "$GIS_OPT_PROJECTION" | tr '[:lower:]' '[:upper:]'`
+   
+   case $PROJECTION in
+      "TRANSVERSE MERCATOR")
+        KNQ_PC=TC
+	KNP_PP="$GIS_OPT_LON_0"
+        KNQ_P1="$KNP_PP"
+	KNQ_P2=1.0        # +k_0  # ??? guess
+	KNQ_P3="$GIS_OPT_LAT_0"  # ??? guess
+	KNQ_P4=NOT_APPLICABLE
+	;;
+      "MERCATOR")
+        KNQ_PC=MC
+	KNP_PP="$GIS_OPT_LAT_0"
+	;;
+      "LAMBERT CONFORMAL CONIC")
+        KNQ_PC=UNKNOWN
+        KNP_PP=UNKNOWN
+	;;
+      "POLYCONIC")
+        KNQ_PC=PH
+	KNP_PP="$GIS_OPT_LON_0"
+        ;;
+      *)
+	g.message -e "Only merc, tmerc, lcc, and poly projections are supported."
+	exit 1
+	;;
+   esac
+
+   DATUM="$GIS_OPT_DATUM"
+else
+   DATUM=`g.proj -p | grep '^datum' | cut -f2 -d: | awk '{print $1}' | tr '[:lower:]' '[:upper:]'`
+
+   PROJ4=`g.proj -j | grep '^+proj=' | cut -f2 -d=`
+   case $PROJ4 in
+      tmerc)
+	PROJECTION="TRANSVERSE MERCATOR"
+	KNP_PP=`g.proj -j | grep '^+lon_0=' | cut -f2 -d=`
+	KNQ_PC=TC
+        KNQ_P1=$KNP_PP
+	KNQ_P2=`g.proj -j | grep '^+k_0=' | cut -f2 -d=`  # ??? guess
+	KNQ_P3=`g.proj -j | grep '^+lat_0=' | cut -f2 -d=`  # ??? guess
+	KNQ_P4=NOT_APPLICABLE
+	;;
+      merc)
+	PROJECTION="MERCATOR"
+	KNP_PP=`g.proj -j | grep '^+lat_ts=' | cut -f2 -d=`
+	KNQ_PC=MC
+	KNQ_P1=UNKNOWN
+	KNQ_P2=$KNP_PP
+	KNQ_P3=NOT_APPLICABLE
+	KNQ_P4=NOT_APPLICABLE
+	;;
+      lcc)
+	PROJECTION="LAMBERT CONFORMAL CONIC"
+	KNP_PP=UNKNOWN
+	#KNQ_P1=UNKNOWN    ???
+	#KNQ_P2=
+	#KNQ_P3=
+	#KNQ_P4=NOT_APPLICABLE
+	;;
+      poly)
+	PROJECTION="POLYCONIC"
+	KNP_PP=`g.proj -j | grep '^+lon_0=' | cut -f2 -d=`
+	KNQ_PC=PH
+	KNQ_P1=$KNP_PP
+	KNQ_P2=UNKNOWN
+	KNQ_P3=NOT_APPLICABLE
+	KNQ_P4=NOT_APPLICABLE
+	;;
+      *)
+	g.message -e "Only merc, tmerc, lcc, and poly projections are supported."
+	exit 1
+	;;
+   esac
+fi
+
+#### set up rest of terms (as best we can)
+#KNP/
+#PP=0        # +lat+0 Projection Parameter (value depends upon Projection) e.g. 135.0
+PI=UNKNOWN   # more proj info?  e.g. 1.0   (scale factor??)   PIxels....?
+# scale : PI ratio.
+#SC=400000,PI=30.000
+#SC=80000,PI=10.000
+#SC=40000,PI=5.000
+#SC=15000,PI=1.000
+#SC=10000,PI=1.000
+# perhaps = (mean_extent / dpi*0.0254)?
+
+
+SP=          # SP=New Jersey,5000   # state plane??
+SK=0.0       # Skew angle e.g. 0.0  (map rotation when north != up)
+TA=90.0      # ?? text angle ?? tranverse angle??  e.g. 90.0
+
+
+# depth unit (FATHOMS) (DX,DY too??)
+UNITS=`r.info -u "$MAP_NAME" | cut -f2- -d= | tr '[:lower:]' '[:upper:]'`
+if [ -z "$UNITS" ] ; then
+   UNITS=UNKNOWN
+fi
+
+# Sounding datum (MEAN LOWER LOW WATER)  `r.info -d $MAPNAME`
+VDATUM=`r.info -d $MAP_NAME | cut -f2- -d= | tr '[:lower:]' '[:upper:]'`
+  # LAT   (lowest astronomical tide)
+  # MEAN LOWER LOW WATER
+if [ -z "$VDATUM" ] ; then
+   VDATUM=UNKNOWN
+fi
+
+
+
+### write out meta-data
+cat << EOF > "$OUTFILE.hdr"
+! Experimental header - NOT FOR NAVIGATIONAL USE - Use at your own risk!
+VER/2.0
+BSB/NA=$MAP_DESCRIPTION
+    NU=$NU,RA=$cols,$rows,DU=$DPI
+KNP/SC=$MAP_SCALE,GD=$DATUM,PR=$PROJECTION
+    PP=$KNP_PP,PI=$PI,SP=$SP,SK=$SK
+    UN=$UNITS,SD=$VDATUM,DX=$ewres,DY=$nsres
+EOF
+
+
+if [ "$PROJECTION" = "TRANSVERSE MERCATOR" ] ; then
+   cat << EOF >> "$OUTFILE.hdr"
+KNQ/SC=$VDATUM,,RM=INVERSE
+    PC=$KNQ_PC,P1=$KNQ_P1,P2=$KNQ_P2,P3=$KNQ_P3,P4=$KNQ_P4
+EOF
+fi
+
+
+### write GCP refs
+if [ "$PROJ_TYPE" -eq 0 ] ; then
+   eval `g.findfile element=group/$GROUP file=POINTS`  #mapset="$GRP_MAPSET"
+   g.region rast="$MAP_NAME"
+   eval `g.region -g`
+   grep -v '^#\|0$' "$file" | \
+     awk -v ROWS=$rows  'BEGIN { N=0; }
+       { N++;
+         printf("REF/%d,%d,%d,%s,%s\n", N, $1, ROWS-$2, $4, $3)}' \
+     >> "$OUTFILE.hdr"
+
+else
+   # use four corners, NW then clockwise
+   eval `g.region -lg`
+   GCP1=0,0,$nw_lat,$nw_long
+   GCP2=$cols,0,$ne_lat,$ne_long
+   GCP3=$cols,$rows,$se_lat,$se_long
+   GCP4=0,$rows,$sw_lat,$sw_long
+   echo "REF/1,$GCP1" >> "$OUTFILE.hdr"
+   echo "REF/2,$GCP2" >> "$OUTFILE.hdr"
+   echo "REF/3,$GCP3" >> "$OUTFILE.hdr"
+   echo "REF/4,$GCP4" >> "$OUTFILE.hdr"
+fi
+
+
+#### write bounding polygon
+if [ -n "$GIS_OPT_BOUNDING_POLY" ] ; then
+   #PLY/1,-36.8899689273687,174.341631581421
+   #PLY/2,-37.1650562040491,174.341510292845
+   #PLY/3,-37.1650709130473,174.763325319593
+
+   POLY_BND="tmp_${GIS_OPT_BOUNDING_POLY}_$$"
+   v.build.polylines in="$GIS_OPT_BOUNDING_POLY" out="$POLY_BND" cats=no --quiet
+   BDY_DEF=`v.out.ascii in=$POLY_BND format=standard | grep -n '^B' | head -n 1`
+   FIRST_LINE=`echo "$BDY_DEF" | cut -f1 -d:`
+   VERTI=`echo "$BDY_DEF" | awk '{print $2}'`
+
+   if [ "$PROJ_TYPE" -eq 0 ] ; then
+      #figure out how to reproject extra x,y points somehow
+     NUM_GCP=`grep -c '^REF/' "$OUTFILE.hdr"`
+     if [ "$NUM_GCP" -ge 6 ] ; then
+        ORDER=2
+     elif [ "$NUM_GCP" -ge 3 ] ; then
+        ORDER=1
+     else
+	g.message -e "Not enough ground control points"
+        exit 1
+     fi
+
+      if [ `echo "$GRASS_VERSION" | cut -f1,2 -d.` = "6.4" ] ; then
+         g.message -e "Simple XY bounding boxes need GRASS 6.5 or newer to transform"
+	 #exit 1
+      fi
+
+      if [ "$MAP_SCALE" -gt 1000000 ] ; then
+	 g.message -w "Double check bounding polygon came out right."
+      fi
+
+      v.out.ascii in="$POLY_BND" format=standard | \
+           sed -e "1,${FIRST_LINE}d" | head -n "$VERTI" | \
+           g.transform group="$GROUP" order="$ORDER" coords=- format="" | \
+           awk 'BEGIN { ROW = 1; } {
+	        if ($1 <= 180.0) {CLON=$1} else {CLON=$1 - 360.0};
+                print "PLY/" ROW "," $2 "," CLON;
+                ROW++;}' >> "$OUTFILE.hdr"
+   else
+      # known map projection
+      v.out.ascii in="$POLY_BND" format=standard | \
+         sed -e "1,${FIRST_LINE}d" | head -n "$VERTI" | \
+         m.proj -od --quiet | \
+         awk 'BEGIN { ROW = 1; } {
+              print "PLY/" ROW "," $2 "," $1;
+	      ROW++;}' >> "$OUTFILE.hdr"
+   fi
+
+   g.remove vect="$POLY_BND" --quiet
+else
+   ### PLY is just region bounding box converted to lat/lon
+   if [ "$PROJ_TYPE" -eq 0 ] ; then
+      if [ `echo "$GRASS_VERSION" | cut -f1,2 -d.` = "6.4" ] ; then
+         g.message -e "Simple XY bounding boxes need GRASS 6.5 or newer to transform"
+	 #exit 1
+      fi
+
+      if [ "$MAP_SCALE" -gt 1000000 ] ; then
+	 g.message -w "Double check bounding polygon came out right."
+      fi
+
+      eval `g.region -g`
+      RESULT=`echo "$w $n|$e $n|$e $s|$w $s" | tr '|' '\n' \
+          | g.transform group="$GROUP" coords=- order="$ORDER" format=""`
+
+      nw_lat=`echo $RESULT | cut -f1 -d' '`
+      nw_long=`echo $RESULT | cut -f2 -d' ' | awk '{if ($1 <= 180.0) {print} else {print $1 - 360.0} }'`
+      ne_lat=`echo $RESULT | cut -f3 -d' '`
+      ne_long=`echo $RESULT | cut -f4 -d' ' | awk '{if ($1 <= 180.0) {print} else {print $1 - 360.0} }'`
+      se_lat=`echo $RESULT | cut -f5 -d' '`
+      se_long=`echo $RESULT | cut -f6 -d' ' | awk '{if ($1 <= 180.0) {print} else {print $1 - 360.0} }'`
+      sw_lat=`echo $RESULT | cut -f7 -d' '`
+      sw_long=`echo $RESULT | cut -f8 -d' ' | awk '{if ($1 <= 180.0) {print} else {print $1 - 360.0} }'`
+      CNR1="0,0,$nw_lat,$nw_long"
+      CNR2="$cols,0,$ne_lat,$ne_long"
+      CNR3="$cols,$rows,$se_lat,$se_long"
+      CNR4="0,$rows,$sw_lat,$sw_long"
+   else
+       # use four corners, NW then clockwise
+      eval `g.region -lg`
+   fi
+
+   echo "PLY/1,$CNR1" >> "$OUTFILE.hdr"
+   echo "PLY/2,$CNR2" >> "$OUTFILE.hdr"
+   echo "PLY/3,$CNR3" >> "$OUTFILE.hdr"
+   echo "PLY/4,$CNR4" >> "$OUTFILE.hdr"
+fi
+
+
+# Set the Map Crosses 180 flag?
+# mc2bsbh/mc2bsbh.cpp figures this out automatically with:
+#   if((maxlon*minlon)<0.0 && (maxlonx < minlonx))
+if [ $GIS_FLAG_D -eq 1 ] ; then
+   echo "CPH/180.0" >> "$OUTFILE.hdr"
+else
+   echo "CPH/0.0" >> "$OUTFILE.hdr"
+fi
+
+
+# datum shift, in *seconds* (floating point), to apply to REF/ points
+# to bring them into line with WGS84.
+# see http://trac.osgeo.org/gdal/ticket/3172
+#   DTM/<latitude shift>,<longitude shift>
+echo "DTM/0,0" >> "$OUTFILE.hdr"
+
+# !#@$!#%!
+if [ -x `which unix2dos` ] ; then
+   unix2dos "$OUTFILE.hdr"
+else
+   if [ `tail "$OUTFILE.hdr" | file - | grep -c CRLF` -eq 0 ] ; then
+      \mv "$OUTFILE.hdr" "$OUTFILE.hdr.unix"
+      #awk 'sub("$", "\r")' 
+      sed -e 's/$/\r/' < "$OUTFILE.hdr.unix" > "$OUTFILE.hdr"
+      \rm "$OUTFILE.hdr.unix"
+   fi
+fi
+
+
+if [ "$GIS_FLAG_E" -eq 1 ] ; then
+   exit 0
+fi
+
+
+if [ "$GIS_FLAG_B" -eq 1 ] ; then
+   r.out.tiff -p in="$MAP_NAME" out="$OUTFILE.tif" compression=packbit
+else
+   if [ ! -x `which convert` ] ; then
+      g.message -e "ImageMagick tools required for raster export. Please install."
+      exit 1
+   fi
+
+   # mmmph bug: ppmquant reduces to 56 colors when I ask for 127
+   # mmmph bug: pnmquant doesn't read from stdin.
+   # mmmph bug: pnmquant uses a ridiculous amount of RAM.
+   # +___________________________________
+   #  use ImageMagick's `convert` instead
+
+   #MAGICK_MEMORY_LIMIT=1500mb
+   #MAGICK_MAP_LIMIT=1500mb
+   #export MAGICK_MEMORY_LIMIT MAGICK_MAP_LIMIT
+
+   r.out.ppm in="$MAP_NAME" output=- | \
+      convert - -colors 127 -compress RLE "$OUTFILE.tif"
+fi
+
+
+if [ "$GIS_FLAG_Z" -eq 1 ] ; then
+   g.remove rast="$GIS_OPT_INPUT.cropped_$$" --quiet
+fi
+
+
+
+g.message "`basename $0` complete."
+exit
+
+# libbsb usage:
+#tif2bsb infile.hdr infile.tif outfile_1.kap
+


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



More information about the grass-commit mailing list