Geostatistics Workshop

FIA ArcView Scripts and Extensions for Geostatistical Analysis

These arcview scripts are all available from various sources, including the ESRI arcscripts page.  Their use can be seen in 2 videos:  loadext.avi and addcoor.avi .  If you already know how to use scripts, you don't need to worry about the videos.  One thing to remember about the scripts:  you need to figure out what needs to be active before running them.  In other words, you need to have your script in one window, and either a view, theme, or table "active" in the other window.  You need to click on this "active" object before running the script.  With extensions, you just need to remember that they belong in the root installation directory wherever you've got arcview installed.  See your system person if you're using UNIX arcview; the scripts belong in the place where you've got the other .avx files stored (in the $AVHOME/ext directory), or, on PC, in the "ext32" directory.

**Note!    You can see how to load and use these scripts via the vidoes "loadext.avi" and "addcoo.avi" on the videos page.

Download Zip File of all Scripts/Extensions Here.

projector.avx:  This extension will ask you for the input and output projections of your point shapefile.  You would use this after 'adding an event theme', using X and Y coordinates in decimal degrees, for example, from the Eastwide Database.  Then, you would project into something like UTM, State Plane, or Albers, to get the points in their correct locations for modelling (there will be some distortion with all of these projections; UTM and State Plane probably have the least).  If you use Arcview 3.2, the projection utility takes the place of this script. ***NOTE:  if you want to change the datum of your data, you need to use an extension like NADCON (from the arcscripts page), or something else (unless using Arcview 3.2).

addcoo.ave:  This script will add x,y coordinates to a point shapefile's attribute table after the shapefile has been projected into a projection of interest (see above).  The new columns in the attribute table will be called "x-coor" and "y-coor".  The arcinfo version of this is the project command with the file option, or the aml pointproject.aml (included with the amls).

xtools.avx:   This extension will offer many arc-info features.  You can also turn a polygon you've created with the drawing tool into a polygon shapefile, and vice versa.  It's a very useful extension, one I highly recommend using.

tools_dd.avx:  This extension will do a lot of things xtools will do, but it will also summarize points inside of polygons, something useful for evaluating the summary statistics from a simulation/kriging run, or the original data.  This script will also create buffers for you!

ungen.ave:  This script will convert a shapefile into "Arc-Info ungenerate format", or a pair of coordinates for each node in a polygon coverage.  The resulting text file is then brought into Surfer, changed slightly, and turned into a Surfer blanking file, or file which can be used to clip out areas of a final map which you don't want to display (e.g., outside of state boundaries.  Very useful.

legendforimage.apr:  This project needs to be opened, and then you can load an arcinfo grid as an image.  Before doing this, you must have run "gridviewer.aml" from the amls below (i.e., you must have a .rmp and .clr file associated with the grid in the same directory as the grid).  The painters palette icon will create a fake shapfile based on that colormap file, and then delete it, leaving the legend, which you can bring into a layout, and thus get a "real" legend based on the colormap of the grid.  This is a way to trick arcview into producing a legend for a grid loaded as an image.

FIA AML's for use with GSLIB and Surfer Geostatistics Programs

These amls are all used for the manipulation of gslib-format output files into Arc-Info and Surfer format grids, or Arcview shapefiles.  They all are cobbled together from my flimsy aml programming background, so they might crash on you mysteriously, and some might fly into never ending loops if you do not enter a parameter correctly, or if you have done a step wrong.  They have worked on all of my datasets, however.

They will interactively ask you for the grid definition, and any other information.  You should have all of these aml's in the same directory, wherever you're doing the work, or you can use the "&amlpath" command in arc to define the amlpath to the directory where you store your aml's (that's what I do; it keeps things neat).

If you do get into an error that keeps recurring, try leaving arc and then coming back in; I don't know the code yet to delete previously defined variables, etc.

***If you get into a horrible never ending loop (you'll know if you do), hit "ctrl-z or ctrl-c" to get out of it.  Then, leave arc, and come back in, and delete any files which were generated.

Download Zip File of all Aml's Here.

postmegridpoint.aml:  Will take your postme output files (17 columns of numbers, with grid specifications the same as gslib output files) and convert them to exported grids (.e00 files), and point shapefiles.  The point shapefiles (which incidentally represent the nodes of the postme output grids) can be loaded into arcview, and the grids can be imported by doing at the arc prompt "import grid .  You can then run "gridviewer.aml" on them to make a raster coverage which you can see in arcview.  Follow the video instructions for this.

**The aml assumes that your output files are named exactly like this:  e m v cv iqr 17 20 30 25 75 83 60 70 80 90 40 10.  If they're not, you need to rename them.

The program will ask you interactively to define your grid.

**Warning!  If it goes crazy, be sure to hit ctrl-z, and then quit out of arc.  I had to eliminate the error handling, or it wouldn't work properly.  So if there's something that's incorrect (your grid definition, missing files, etc.), it will crash horrifically!**

postmegrid.aml:  This is a good one in that it asks you for a boundary file, and the grid definition, and it will summarize the zones defined by your boundary file (for example, counties), and give you a statistics file.  For example, you can calculate the mean and standard deviation of estimates inside county "A" from each postme output percentile map, as well as the mean and sd of the original data inside of county A (in arcview or excel), and compare the means and sd's, and decide how close your estimates are to your original data (if that's your goal). You might want to minimize the difference between your county level estimates of mean and sd and the original data values of mean and sd.  You can choose a given postme output map, say the 60th percentile, if its mean and sd is most similar to those of all of the original data, across all counties.  See the videos for this.

The use of this .aml is complicated by the fact that the "zone" coverage must be a true arcinfo coverage (not a shapefile), and have a field which is unique and which has an integer value.  (You can turn your zone shapefile into a coverage by running the shpoly.aml, included with the aml's.)  You can check to see if your zone id is integer by doing "items <zonecov>.pat", and see if your FIPS (census county id) code, for example, is of type I.  If it's not, you need to go into INFO and "alter" the item type for your zone id, such as fips, to "I".  If you don't know how to do this, ask a GIS person or feel free to call me.  See the "fipser.aml" program to see the steps needed for doing this.

Your output from this will be a bunch of zipped up grids and a file called "allstats.txt".  This file needs to be opened and altered in Excel before completing the statistics comparison.  Definitely see the video about this.  Reading the aml will be helpful, too.  The aml takes the column of numbers, adds a header which it creates by prompting you for the grid definition, turns this into a grid, then turns the grid upside down so that it's correctly georeferenced, then it turns the polygon zone coverage into a grid, and summarizes your postme grid based on the "zone" grid you just created.  This zone grid is just a grid file with each pixel falling in a given county, e.g., having the id of that county as its value.  See the arcinfo help for zonalstats and polygrid for more details on this.

gridder.aml:  Turns a single headerless gslib file (just one--a file with a column of numbers, each one representing a grid node) into both an arcinfo grid and a Surfer file grid  (.grd -- ascii grid with header).  It prompts you interactively for the grid definition.

stacker.aml:  This is an unfinished aml that will take an sgsim output file and turn it into a stack of grids (if I ever finish it).  I'm providing it here so you can get some ideas about the use of the "sed" command in unix.  (you can use it to delete the header off a text file).  The concept is:  1. attach a header to the large column of stacked simulations (from sgsim); 2. then, do asciigrid, specifying the extent of just one iteration (1 layer), flip this grid; 3. Now, delete away the top nrows*ncols + 6 rows from the file using the sed unix command, and repeat n times, n being the number of realizations.  Then, you'll have 1 grid for each realization, and make comparisons of multiple realizations.  This works, but for some dumb reason I deleted the meat of the aml.  : (

postmesurfer.aml:  This will take a postme output file and turn it into just a Surfer Grid.  You need to edit the aml to define the z lo and the z hi (see the Surfer help for ascii grid definition), depending on your data.

gridviewer.aml:  Will prompt you for an arcinfo grid, ask you about the legend definition (colors and levels of the legend),  show you the classified grid, and write a colormap (.clr) file and a remap table (.rmp) based on your entries.  These can be used in arcinfo to make a real legend, or the .clr file can be read by arcview if you load the grid as an image.  To properly view the colors in arcview, you need to double click the grid_loaded_as_image in the table of contents, and choose "interval", and choose the proper # of levels you entered in gridviewer. as the # of intervals.  Then, to make a legend, you can either create your own legend, or create a legend using the "legendforimage.apr procedure", where you make a fake shapefile with a legend, and then bring the fake legend into a layout, edit the legend, and then you're done.  To do this, you must open your grid-image from inside this project, which is included.  See the videos for details.

To avoid using Arcview for projection:

pointproject.aml:  This will take a text file of column of long, then lat, in decimal degrees, separated by spaces, and project it into Albers NAD27 meters.  It assumes that your long and lat are in decimal degrees, Nad27, e.g., -76.3343 43.2222.  Note:  you will need a file called geoalb.prj in the same directory as the aml and the text file for this to work.  The file is included with the aml's.  Your file should be like this:
long lat
long lat
etc.
You will get a file like this:
albx alby
albx alby
etc.

The Albers projection parameters can be seen, and can be seen in the second half of the "geoalb.prj" file, included.

textpoint.aml:  This will take a "generate format" text file and turn it into a shapefile. It allows you to avoid using arcview to make a point shapefile.  A generate format text file is like this:
unique_id1,x,y
unique_id2,x,y
unique_id3,x,y
.......
unique_idn,x,y
END

This step should be done with the x and y from the results of pointproject.aml.  It's really just as easy to do this in Arcview, however!

Previous: Videos