[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