Post by Malte ThomaDoes a conversion of the BEDMAP
http://www.antarctica.ac.uk/aedc/bedmap/
datasets to the GMT format exists?
Greetings,
Malte
Hi Malte,
Try the bourne shell script below. This simply makes a grdfile out of the
projected data. To plot the data with lon/lat gridlines overlaid, you
would do something along the lines of:
1. plot the projected data with -JX$width
2. overplot the gridlines (or any other data you have in lon/lat form)
with -JS0/-90/$width, with a region that matches the projected grid
exactly.
(By the way, I generally use perl for non-trivial GMT scripts these days,
since all the calls to grep, cut, bc, etc. get kinda ugly, not to mention
slow.)
Bruce
--
Bruce Raup Phone: 303-492-8814
National Snow and Ice Data Center, U. of Colorado, 449 UCB, Boulder, CO 80309
#!/bin/sh
#
# $Id: asc2grd,v 1.4 2001-08-24 15:54:12-06 braup Exp braup $
#
# To create GMT grdfiles from the BEDMAP ASCII files:
if [ $# != 1 ]; then
echo "Usage: $0 infile"
echo "where infile is a ASCII file from the BAS BEDMAP project"
exit 1
fi
infile=$1 # edit this as you please
if [ ! -e $infile ]; then
echo "File $infile not found"
exit 1
fi
grdfile=`basename $infile .asc`.grd
psfile=`basename $infile .asc`.ps
cptfile=`basename $infile .asc`.cpt
echo "Input file: $infile"
echo "Grid file: $grdfile"
echo "PS file: $psfile"
num_hdr_lines=6
head -$num_hdr_lines $infile > hdr_$$
tail +$num_hdr_llines $infile | tr -s ' ' '\n' > tail_z_$$
cat hdr_$$ tail_z_$$ > ${infile}.zcol
rm tail_z_$$
nc=`grep ncols hdr_$$ | tr -s ' ' ' ' | cut -d" " -f2`
nr=`grep nrows hdr_$$ | tr -s ' ' ' ' | cut -d" " -f2`
cellsize=`grep cellsize hdr_$$ | tr -s ' ' ' ' | cut -d" " -f2`
xll=`grep xllcorner hdr_$$ | tr -s ' ' ' ' | cut -d" " -f2`
yll=`grep yllcorner hdr_$$ | tr -s ' ' ' ' | cut -d" " -f2`
nodata=`grep NODATA_value hdr_$$ | tr -s ' ' ' ' | cut -d" " -f2`
rm hdr_$$
echo "Number of columns: $nc"
echo "Number of rows: $nr"
echo "cell size: $cellsize"
echo "xllcorner: $xll"
echo "yllcorner: $yll"
echo "NODATA value: $nodata (will be converted to NaNs in $grdfile)"
echo "Creating $grdfile ..."
# create region
W=$xll
S=$yll
E=`echo "scale=10; $xll + ($nc-1)*$cellsize" | bc`
N=`echo "scale=10; $yll + ($nr-1)*$cellsize" | bc`
region="$W/$E/$S/$N"
xyz2grd ${infile}.zcol -R$region -I$cellsize -ZTLa -N$nodata -G$grdfile -V
echo "Done"
minmax ${infile}.zcol
rm ${infile}.zcol
echo "Creating PostScript display of $grdfile ..."
#grd2cpt $grdfile -Cglobe -S0000/4000/100 -Z -V > $cptfile
grd2cpt $grdfile -Cglobe -Z -V > $cptfile
#grd2cpt $grdfile -Cglobe -Z -V > $cptfile
xdim=17 # cm
ydim=`echo "scale=5; $xdim*$nr/$nc" | bc`
grdimage $grdfile -C$cptfile -JX${xdim}c/${ydim}c -R$region -B500000 -V > $psfile
echo "Done"
# End of script
