Add a column to a text file from a raster

Text files are great

I have mixed feelings about software interoperability. Sure GRASS interfaces with R, and R can interface with Excel and Access if you have the ODBC drivers. But, often users do not have the energy to learn the best techniques for moving data back and forth between programs. Further, there always seems to be some quirk with data that prevents the transfer from going as smoothly as the user's manual examples do. Then, I remember the lowly text file. Whether you use comma separated values (csv) files or tab-delimited files, every software platform on the planet can read text files of records (lines) and fields (values on a line separated by commas,tabs or some other character). If you're trying to send data to someone else, you can be confident that your recipient can read them. It may not be the flashiest or fastest way to move data between people and systems, but it always works. Finally, when you use plain text files for data, you have many more options for working with the files. Tools such as bash, awk, sed, perl, and any other programming language you can think of will can be brought to bear on data analysis problems.

Attaching raster data to points in GIS

Switch gears for a moment from a common file format to a common analysis task. For researchers who work with spatial point patterns (like plots), a regular step in data analysis is to pull data from raster layers where they intersect the point data and add it as columns to the point data. Informally, these are known as "drill-down" operations. You have a stack of maps and you want to drill down through them vertically at point locations, grab the data, and add it to the points. This is an a relatively simple task in ArcGIS if you use Hawth's tools. It's also pretty easy in GRASS, though it involves a couple of steps. The function v.what.rast will put raster values into a column of a point dataset. The column must already exist in the dataset, which can be done through the command v.db.addcol and v.what.rast function just fills in the missing values with values where the raster intersects the points.

Adding raster values to text files

Now, back to text files. I recently worked with a graduate student who wanted elevation values from a lidar-derived digital elevation model (dem) on 4 samples totalling about 1100 plots. She is not a GRASS user, and so she sent me the locations of all the plots along with an identifier in csv format. I wanted to send her back the same csv files with elevation simply added as a field. I wrote a bash script to accomplish this quickly and repeatably for her files. In this script, all the GRASS functions are hidden from the user. Users simply input a csv file and the raster column you wish to add, and a csv file is output with that column added. The script requires the awk script that I posted previously, which recognizes the types of varibles present in the input csv. I have included a portion of the full script below. GRASS has a very particaular code style that it uses for running scripts natively (See g.parser for details). This makes the script very long. So, I've cut out the not-so-important parts. If you would like the full file you can get it here.

v.ascii.addrastcol

#!/bin/sh
##################################################################
#
# MODULE:       v.ascii.addrastcol
# AUTHOR(S):    R. Todd Jobe 
# PURPOSE:      Add a raster value field to a csv of points
# COPYRIGHT:    (C) 2009 by R. Todd Jobe
#
#  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.
#
##################################################################
# . . . Grass header information removed
# . . . The cut code creates variables GIS_OPT_* and GIS_FLAG_*
# . . . Get the full code at:
# . . . http://www.unc.edu/~toddjobe/src/v.ascii.addrastcol
##################################################################
# Need space-delimited arguments for looping
GIS_OPT_INPUT=${GIS_OPT_INPUT//,/ }

# test for input raster map
eval `g.findfile element=cell file="$GIS_OPT_RAST"`
if [ ! "$file" ] ; then
   g.message -e "Raster map <$GIS_OPT_RAST> not found in mapset search path"
   exit 1
fi

for i in $GIS_OPT_INPUT; do
   # Check that input files exist
   if [ ! "$i" ] ; then
 echo "file $i does not exist" >&2
 exit $E_NOTFOUND
   fi

   # Create header line if necessary and get variable types
   if [ $GIS_FLAG_H -eq 1 ]; then
 skip=1
   else
 skip=0
   fi

   # Intermediate variables
   map=`basename $i .csv`
   tmp=`tempfile -s ".csv"`
   of="${map}${GIS_OPT_SUFFIX}"
   col=`typeVars.awk -v h=$GIS_FLAG_H -F, $i`

   # Get the raster type
   if `r.info -t $GIS_OPT_RAST | grep -q "FCELL"`; then
 rtype="double precision"
   else
 rtype="int"
   fi
  
   # Grass-specific calls
   v.in.ascii --overwrite input="$i" output="$map" fs="$GIS_OPT_FS" \
 skip=$skip column="$col" x="$GIS_OPT_X" y="$GIS_OPT_Y"
   v.db.addcol map="$map" \
 columns="$GIS_OPT_RAST $rtype"
   v.what.rast vector="$map" column="$GIS_OPT_RAST" \
 raster="$GIS_OPT_RAST"
   v.out.ascii input="$map" fs="," \
 columns="$GIS_OPT_RAST" | \
   awk -F, -v meters="$GIS_FLAG_M" -v head="$GIS_FLAG_H" \
 -v rast="$GIS_OPT_RAST" '
   BEGIN{ if(head == 1) print rast;}
   {
        # Convert Ft. elevation to meters if necessary
        if($4!="" && meters==1) $4=$4 * 0.30480060
        print $4
   }' >"$tmp"
   # Join raster values to input and write to output
   paste -d$GIS_OPT_FS "$i" "$tmp" >"$of"

   # Cleanup
   rm "$tmp"
   if [ $GIS_FLAG_R -eq 1 ] ; then
 g.remove -f --quiet vect=$map
   fi
done

To see the useage for the full script, you can simply start GRASS and from the command line type v.ascii.addrastcol --help. Here's a short example for a random set of 25 points on Ft. Bragg North Carolina.

GRASS 6.4.0RC3 (ftBragg):~ > head locs.csv
"","V1","V2","V3"
"1",594927.510605283,153893.107849400,1
"2",590884.780751329,155855.275043186,2
"3",594095.580107841,155284.234116311,3
"4",605646.865006463,154397.274762958,4
"5",603002.815879479,153284.801461238,5
"6",603453.085431226,152633.248500899,6
"7",602981.046299442,152156.288326646,7
"8",604160.50828444,151921.617483228,8
"9",603752.189821169,152062.642227756,9
GRASS 6.4.0RC3 (ftBragg):~ > v.ascii.addrastcol --quiet -hr \
input=locs.csv rast=iLandcover x=2 y=3 --quiet -r
GRASS 6.4.0RC3 (ftBragg):~ > head locs_isect.csv
"","V1","V2","V3",iLandcover
"1",594927.510605283,153893.107849400,1,3
"2",590884.780751329,155855.275043186,2,3
"3",594095.580107841,155284.234116311,3,1
"4",605646.865006463,154397.274762958,4,3
"5",603002.815879479,153284.801461238,5,3
"6",603453.085431226,152633.248500899,6,3
"7",602981.046299442,152156.288326646,7,3
"8",604160.50828444,151921.617483228,8,1
"9",603752.189821169,152062.642227756,9,3