[GRASSLIST:6624] Grass GIS and Fragstats on Linux!

[Sander Oom]

Dear Fragstats and Grass GIS users,

To my pleasant surprise I have successfully run Fragstats on a Linux 
machine in conjunction with my Grass GIS database.

I tried the landscape index tools r.le within Grass, but ran into 
several problems (see earlier posts on the Grass mailing list). Than I 
started looking for an alternative. Fragstats is a powerful and well 
tested program, but it is only available for Windows.

However, when running the Fragstats setup program using CrossOver Office 
on my Linux machine, the installer finished without any errors! I have 
tried several options and methods and have not seen the program crash 
yet! I guess the success is a result of the combination of the little 
dependency of Fragstats on odd windows dll's and the progress made in 
the Wine project, the basis for CrossOver Office. I have not tested 
Fragstats on Wine.

Below you find the Grass code to prepare your data for Fragstats. I used 
the Albers Equal Area projection to prepare my rasters. Using sed (Linux 
utility) I removed the header from the exported ascii file. Fragstats 
only wants the points. I put the header content in a file called 
'asc_header_aea.txt' so I could paste it back in later, as Fragstats 
returns exactly the same grid of points. Fragstats creates a new folder 
with an ascii grid (no header nor an extension!) for each index 
calculated, so I wrote a loop to go through the files, add the header 
and read the result into Grass. Then I queried the raster with a points 
vector and exported the results to a text file.

Hope this is helpful for other Grass/Linux users!

Cheers,

Sander.

########################################################################
## Export to ASCII and prepare file for Fragstats
dirOut=/grassdata/output
rastIn="nlc31"
rastOut=$dirOut/nlc31_aea
r.out.arc input=$rastIn output=$rastOut.asc dp=6
## remove header and replace -9999 with 9999
cd $dirOut
sed '1,6d' < nlc31_aea.asc > nlc31_aea_frgs.asc

########################################################################
## Import Fragstats results
## query with sa10min
## export to text file for analysis in R

map=nlc31
   for index in MSIDI MSIEI SHDI SHEI SIDI SIEI ; do

     ## set resolution to 1 km; projection in meters!
     g.region res=1000

     ## Add header to fragstats result and paste into ASCII file
     frgsIn=/grassdata/output/${map:-nlc31}_aea_frgs.asc_MW1/${index:-MSIDI}
     hdrAEA=/grassdata/output/asc_header_aea.txt
     ascOut=/grassdata/output/${map:-nlc31}_aea_${index:-MSIDI}.asc
     echo $frgsIn
     cat $hdrAEA $frgsIn > $ascOut

     ## Import ASCII
     ascIn=$ascOut
     rastOut=${map:-nlc31}${index:-MSIDI}
     r.in.arc input=$ascIn output=$rastOut

     ## Query raster with sa10min and export to file
     vectTmp=tmp
     g.remove vect=$vectTmp
     g.copy vect=sa10min,$vectTmp
     dirDBF="/grassdata/sa-aea/macwits/dbf/"
     db.columns table=$vectTmp driver=dbf database=$dirDBF
     ## Add new column to attribute table
     echo "ALTER TABLE $vectTmp ADD COLUMN value double " | db.execute
     db.columns table=$vectTmp driver=dbf database=$dirDBF
     ## Query raster with vector
     v.what.rast vect=$vectTmp rast=$rastOut col=value
     dirFrags=/data1/gis/fragstats/macwits/output
     textOut=$dirFrags/${map:-nlc31}_sa10min_${index:-MSIDI}.txt
     echo " SELECT * FROM $vectTmp " | db.select fs="," > $textOut

   done
########################################################################


-- 
--------------------------------------------
Dr. Sander P. Oom
Animal, Plant and Environmental Sciences,
University of the Witwatersrand
Private Bag 3, Wits 2050, South Africa
Tel (work)      +27 (0)11 717 64 04
Tel (home)      +27 (0)18 297 44 51
Fax             +27 (0)18 299 24 64
Email   sander at oomvanlieshout.net
Web     www.oomvanlieshout.net/sander
---------------------------------------------