Abstract:
A collection of articles, methods, and examples of how open source software tools can assist with research in the field of soil science. Additional links to methods in environmental sciences within the realm of statistics, interpolation, morphometrics, and more.
UCD Resources:
Garrett and I are leading a seminar-style course covering some of this material, Fall quarter, 2006.
Articles:
Reference Material
Some initial hints of use of various scripting languages.
Python
Ruby
Premise:
Thought it would be fun to compute how much ink a given poster requires, per unit area of paper, when sending to the department large-format printer. The Python Imaging Library provides several modules suitable for low-level operation on image data. A simple (and probably very inefficient) script was developed to compute the white/black percentage of an image. A script like this could be used to adjust a per-poster "ink cost", which would hopefully prevent people from wasting ink. Obviously, this computation is scale-dependent, so standardized rasterization parameters would have to be set in order for the "ink cost" calculation to be fair. More generalized or efficient approaches are always welcomed.
Implementation: (when copying/pasting note whitespace in blocks)
Premise
Wanted a simpler way to access the USGS seamless elevation look-up service. Python seemed like a logical start. Note that the response from the USGS webservice is not correctly identified as valid XML by the python XML-parser. Therefore there is a small amount of scrubbing used to coerce the response into valid XML. Comments on why this is, or is not, a good idea are welcome. Here is a link to an interactive version.
Update It looks like the USGS service does not accept POST-style requests. I have made some small changes to the script below.
Example Listing
Incantation
Input File CA cities in WGS84 (lon,lat,id) records
-122.32,37.78,Alameda NAS -120.53,41.48,Alturas -124.1,40.98,Arcata -119.05,35.43,Bakersfield -121.45,39.13,Beale AFB -116.95,33.93,Beaumont -116.62,35.28,Bicycle Lk -116.68,34.27,Big Bear Apt -118.6,37.60,Bishop -120.7,39.28,Blue Canyon -114.72,33.62,Blythe -118.37,34.20,Burbank
Output
-122.320000,37.780000,Alameda NAS,2.123548,NED Contiguous U. S. 1/3W arc second elevation data -120.530000,41.480000,Alturas,1331.056396,NED Contiguous U. S. 1/3W arc second elevation data -124.100000,40.980000,Arcata,63.749836,NED Contiguous U. S. 1/3W arc second elevation data -119.050000,35.430000,Bakersfield,148.473618,NED Contiguous U. S. 1/3W arc second elevation data -121.450000,39.130000,Beale AFB,29.950783,NED Contiguous U. S. 1/3W arc second elevation data -116.950000,33.930000,Beaumont,792.840881,NED Contiguous U. S. 1/3W arc second elevation data -116.620000,35.280000,Bicycle Lk,711.765869,NED Contiguous U. S. 1/3W arc second elevation data -116.680000,34.270000,Big Bear Apt,1697.037720,NED Contiguous U. S. 1/3W arc second elevation data -118.600000,37.600000,Bishop,2209.536865,NED Contiguous U. S. 1/3W arc second elevation data -120.700000,39.280000,Blue Canyon,1603.887573,NED Contiguous U. S. 1/3W arc second elevation data -114.720000,33.620000,Blythe,120.626007,NED Contiguous U. S. 1/3W arc second elevation data -118.370000,34.200000,Burbank,224.287033,NED Contiguous U. S. 1/3W arc second elevation data -117.350000,33.300000,Camp Pendlet,21.895588,NED Contiguous U. S. 1/3W arc second elevation data
Working with DBF files
http://www.graphviz.org/
Example DOT Code:
Example DOT Code:
A collection of links to documents which I have found helpful in working with the LaTeX document preparation system. An internal link to all documents related to LaTeX can be accessed here.
Getting Started
CV Creation
Journal Submission
Collaborative Writing Ideas
Premise:
See attached files for a bibliography style (.bst) file that is compatible with the SSSAJ (Soil Science Society of America Journal) style guidelines.
General Notes:
Documentation:
http://www.andy-roberts.net/misc/latex/latextutorial3.html
A Custom Bib Style:
Conversion of a .tex document to open document format (RTF)
Conversion of a .tex document to open document format (ODT)
Some notes:
Conversion of PDF, EPS, or Postscript figures to EMF (enhanced metafile format) for windows apps.
Note that these instructions are for unix-like systems. EMF support can be added to pstoedit via libEMF. The libEMF library contains several windows-specific coding practices, along with some sloppy use of include files. Instructions for compiling with a modern version of GCC (4.x) can be found here.
Premise:
The standard Latex table environment can be difficult to extend, especially when one wants to use modern constructs such as table foot notes, etc. The ctable package is a convenient approach to solving this problem. A complete worked example of a table typset using the \ctable command is presented below. Note that \ctable is a command, and therefore does not allow blank newline characters. One way to maintain readability is to trick Latex by adding a comment character to all blank lines. Link to PDF manual for ctable.
Example of table produced with ctable
Code used to produce the example above. This example was produced with TexLive.
%
% start the table: note that we cannot have extra newline characters in the ctable defs
\ctable[
cap = {logistic regression parameters},
%
caption = {Logistic regression model (M3) parameters. Coefficients, standard error, z-values, and p-values are included for each term used in the model. A separate slope and intercept term was fit to each geologic class.},
%
label = {aspect_effect:table:glm_model_properties},
%
%
]{lccc}{
%
\tnote[$\ast$]{Marginal p-values are used to determine whether each term is significantly different than 0.}
%
}{ \FL
%
Model Term & Value & Std. Error & Marginal p-values\tmark[$\ast$] \ML
%
Intercept & & & \\ \cmidrule(r){1-1}
%
\;\; andesite & 1.146e+01 & 6.920e+00 & 0.09781 \\
\;\; clastic\_volcanic & 3.471e+00 & 2.234e+00 & 0.12024 \\
\;\; coarse\_sedimentary & 4.446e+00 & 2.390e+00 & 0.06285 \\
\;\; felsic\_intrusive & 5.384e+00 & 2.714e+00 & 0.04725 \\
\;\; fine\_sedimentary & 8.527e+00 & 4.013e+00 & 0.03361 \\
\;\; rhyolite & 1.075e+01 & 4.071e+00 & 0.00828 \\
\;\; tuff & 1.657e+01 & 9.116e+03 & 0.99855 \\
%
%
\\
Slope & & & \\ \cmidrule(r){1-1}
%
\;\; andesite & -1.356e-03 & 8.448e-04 & 0.10858 \\
\;\; clastic\_volcanic & -3.588e-04 & 2.851e-04 & 0.20825 \\
\;\; coarse\_sedimentary & -8.165e-04 & 3.336e-04 & 0.01438 \\
\;\; felsic\_intrusive & -6.352e-04 & 3.503e-04 & 0.06977 \\
\;\; fine\_sedimentary & -1.122e-03 & 5.139e-04 & 0.02897 \\
\;\; rhyolite & -1.361e-03 & 4.873e-04 & 0.00522 \\
\;\; tuff & 6.155e-14 & 1.200e+00 & 1.00000
%
%
\LL}
PostGIS is a set of extensions to the relational database management system PostgreSQL, which provide access to spatial constructs, operators, and functions. In other words, PostGIS is a spatially-aware version of Postgresql. Paul Ramsey of Refractions Research has written up a nice summary of how were are using PostGIS.

Tuning Tips:
PostGIS Syntax Examples:
PostGIS Presentations:
GDAL/OGR tools This approach allows simultaneous conversion of coordinate systems, but is less flexible with respect to generation of new tables in PostGIS.
ogr2ogr -f "PostgreSQL" PG:'dbname=ssurgo_combined user=xxxx password=xxxx host=postgis.server.edu' input_file.shp
ogr2ogr output_file.shp PG:'dbname=ssurgo_combined user=xxxx password=xxxx host=postgis.server.edu' tablename
Note that tables must be correctly 'registered' in the geometry_columns table for this to work:
INSERT INTO geometry_columns VALUES ('','public','tablename','wkb_geometry',2,SRID,'geomtype');
PostGIS Loader/Dumper This approach is the simplest, but does not allow on-the-fly conversion of coordinate systems.
shp2pgsql -s SRID -c -g wkb_geometry -I shapefile.shp schema.table | psql -U username -h host database
Note that SRID is the PostGIS 'spatial ref. sys. id' (see the spatial_ref_sys table). See the manual page for shp2pgsql for a complete list of arguments and their meanings.
pgsql2shp -f shapefile.shp -h host -u username -P password -k -g wkb_geometry database schema.table
See the manual page for pgsql2shp for a complete list of arguments and their meanings.
Where tablename is your newly created table, SRID is the SRID (spatial reference ID) for the geometry in this table, and geomtype is the type of geometry: POINT, LINE, POLYGON, etc.
CSV format, from within the psql client
\copy tablename TO 'filename.csv' CSV HEADER
CSV format, via psql client
echo "select column_list from table_list " | psql --tuples --no-align -F "," database > file.csv
Tabular data to HTML format, via psql client See output below:
echo "select column_list from table_list " | psql --html database > file.html
HTML output from psql
| area | compname |
|---|---|
| 132472.230854819 | Hilmar variant |
| 322819.967391312 | Oneil |
| 362729.418301135 | Carranza |
| 431948.171760353 | Tuff rockland |
| 448784.927049035 | Gravel pits |
| 500763.225267798 | Snelling variant |
| 518860.954990617 | Foster |
| 571640.132661382 | Alamo |
| 648973.748756059 | Toomes |
| 924327.631201791 | Dumps |
(10 rows)
Overview
The ST_Affine() function from PostGIS is useful for manipulating geometries, but requires the elements of a transformation matrix. This page documents progress on automating the computation of the transformation matrix by least-squares (Bruce Rindahl) via SQL. This would allow a pure PostGIS solution to computing and applying affine transformations to geometry data.
An open-source algorithm for computing the transformation matrix
Example code from GRASS (v.transform) was used as a template.
Approach
Compute transformation matrix based on a table of control points, stored as numbers.
Evaluation of results
Comparable to output from a similar analysis done in R, and the original algorithm as implemented in v.transform (GRASS).
Example of bad Tiger data in Stanislaus County: Red lines are the original road network, green lines are the corrected road network.
The Problem
The US Census does a nice job of collecting all sorts of geographic and demographic information every 10 years. This data is available free of charge in the rather complex and soon to be replaced TIGER/LINE format. While this data covers the entire US down to the local road level, there are numerous errors and even extreme cases of coordinate-shift. Here is an example from Stanislaus County, California. The original TIGER data (red lines) are offset several hundred meters from the imagery. While it is not clear what may have caused the problem, it can be fixed without much effort using an affine transformation. We do not have the transformation matrix, however it can be 'fit' to a set of control points by several methods. The general form of the affine transform can be conveniently represented in homogeneous coordinates as:
Affine Matrix Form in homogeneous coordinates:New coordinates on the left-hand side, old coordinates on the right-hand side. The transformation matrix is the 3x3 matrix in the center.
The Solution
We first need a set of control points, good and bad coordinates. This can be accomplished in several ways, we used the d.where command in GRASS:
Computing the transformation matrix can be done with a simple regression between 'good' and 'bad' coordinates in R. Note that this approach was suggested by Prof. Brian Ripley on the R-help mailing list.
Compute the Affine Transformation Matrix in R
Establishing the transformation based on control points: Red points represent where the coordinates should be. Black points are the original and incorect coordinates.
Check Affine Transform Results in PostGIS
Perform Affine Transformation in PostGIS
Regression Diagnostics
Response nx :
Call:
lm(formula = nx ~ x + y, data = d)
Residuals:
Min 1Q Median 3Q Max
-207.088 -23.856 8.614 21.245 161.610
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.017e+03 1.369e+03 3.666 0.00079 ***
x 1.002e+00 6.654e-04 1506.386 < 2e-16 ***
y 9.190e-03 9.419e-04 9.756 1.20e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 52.25 on 36 degrees of freedom
Multiple R-Squared: 1, Adjusted R-squared: 1
F-statistic: 1.322e+06 on 2 and 36 DF, p-value: < 2.2e-16
Response ny :
Call:
lm(formula = ny ~ x + y, data = d)
Residuals:
Min 1Q Median 3Q Max
-39.835 -18.459 -4.556 15.311 94.226
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.814e+04 7.409e+02 -37.98 <2e-16 ***
x -1.352e-02 3.602e-04 -37.54 <2e-16 ***
y 9.974e-01 5.099e-04 1956.22 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 28.28 on 36 degrees of freedom
Multiple R-Squared: 1, Adjusted R-squared: 1
F-statistic: 2.187e+06 on 2 and 36 DF, p-value: < 2.2e-16
First the output from R:
Looking at the residuals from the regression model used to map bad coordinates (x,y) to good coordinates (nx,ny):
x y nx ny resid
1 -2078417 -14810.570 -2078314 -14838.378 46.617600
2 -2078743 -16057.955 -2078636 -16081.790 62.041274
3 -2077261 -16435.348 -2077170 -16463.156 40.905132
4 -2076709 -14405.369 -2076606 -14433.177 29.406399
5 -2074179 -15830.248 -2074084 -15901.558 33.111981
6 -2073850 -15707.435 -2073763 -15798.554 37.362736
7 -2073450 -13873.171 -2073359 -13920.712 21.623235
8 -2072359 -15204.613 -2072276 -15323.138 38.997678
9 -2072545 -14402.596 -2072450 -14513.219 32.918889
10 -2072189 -16022.434 -2072098 -16129.106 33.834074
11 -2071991 -16856.058 -2071928 -16942.976 6.277554
12 -2068407 -12999.396 -2068296 -13133.170 6.579285
13 -2069870 -12613.813 -2069764 -12731.848 2.357631
14 -2067635 -13188.253 -2067517 -13337.765 11.604519
15 -2066931 -13377.110 -2066809 -13518.753 22.719625
16 -2067411 -15084.692 -2067313 -15190.924 41.907273
17 -2066795 -18714.093 -2066741 -18846.019 14.541358
18 -2066384 -17080.538 -2066299 -17212.464 26.717495
19 -2068634 -19742.339 -2068580 -19835.464 27.483654
20 -2053326 -16930.710 -2053276 -17226.351 65.746074
21 -2051797 -17321.500 -2051899 -17579.762 227.516944
22 -2068307 2826.921 -2068066 2638.276 12.587853
23 -2067543 2648.205 -2067328 2449.631 37.729747
24 -2067126 4276.510 -2066904 4081.246 46.774630
25 -2066748 4170.604 -2066527 4001.816 59.509843
26 -2066068 2292.295 -2065860 2094.699 46.681553
27 -2065337 2107.872 -2065126 1900.397 43.386956
28 -2064606 1913.570 -2064378 1692.922 26.460389
29 -2064199 3558.561 -2063961 3356.401 47.742696
30 -2037464 6512.455 -2037076 5864.398 50.762994
31 -2036722 6825.682 -2036338 6199.227 22.699467
32 -2036876 6366.642 -2036498 5742.888 22.120176
33 -2040225 7150.180 -2039706 6575.029 161.631199
34 -2041064 7144.779 -2040732 6569.629 26.657582
35 -2044702 -15548.033 -2044564 -16024.903 23.844817
36 -2043992 -15723.521 -2043824 -16223.282 48.063840
37 -2043790 -14907.119 -2043611 -15383.990 34.844851
38 -2040616 -14820.445 -2040453 -15349.974 21.196233
39 -2039595 -15081.427 -2039485 -15588.263 47.263287
The Root-Mean-Square-Error (RMSE) for the fitted transform (in meters) is:
The output from v.transform on the same set of control points:
Transformation Matrix | xoff a b | | yoff d e | ------------------------------------------- 5301.399323 1.002469 0.009172 -28155.882288 -0.013530 0.997547 -------------------------------------------
full output including the residuals:
CHECK MAP RESIDUALS
Current Map New Map
POINT X coord Y coord | X coord Y coord | residuals
1. -2078417.36 -14810.57 | -2078314.07 -14838.38 | 46.81
2. -2078743.11 -16057.95 | -2078635.85 -16081.79 | 62.22
3. -2077261.34 -16435.35 | -2077169.97 -16463.16 | 41.05
4. -2076709.16 -14405.37 | -2076605.87 -14433.18 | 29.59
5. -2074178.76 -15830.25 | -2074083.67 -15901.56 | 33.21
6. -2073849.93 -15707.44 | -2073762.78 -15798.55 | 37.42
7. -2073449.80 -13873.17 | -2073358.68 -13920.71 | 21.62
8. -2072358.86 -15204.61 | -2072275.89 -15323.14 | 39.02
9. -2072544.55 -14402.60 | -2072449.73 -14513.22 | 32.97
10. -2072188.97 -16022.43 | -2072098.11 -16129.11 | 33.87
11. -2071991.43 -16856.06 | -2071928.22 -16942.98 | 6.27
12. -2068406.55 -12999.40 | -2068296.38 -13133.17 | 6.60
13. -2069870.19 -12613.81 | -2069763.96 -12731.85 | 2.33
14. -2067635.38 -13188.25 | -2067517.35 -13337.76 | 11.63
15. -2066931.10 -13377.11 | -2066809.13 -13518.75 | 22.74
16. -2067411.11 -15084.69 | -2067312.75 -15190.92 | 41.93
17. -2066795.16 -18714.09 | -2066740.84 -18846.02 | 14.64
18. -2066383.87 -17080.54 | -2066298.50 -17212.46 | 26.74
19. -2068634.37 -19742.34 | -2068580.05 -19835.46 | 27.53
20. -2053326.48 -16930.71 | -2053275.51 -17226.35 | 66.09
21. -2051797.30 -17321.50 | -2051899.25 -17579.76 | 227.91
22. -2068307.24 2826.92 | -2068065.64 2638.28 | 12.41
23. -2067542.73 2648.21 | -2067327.61 2449.63 | 37.44
24. -2067125.72 4276.51 | -2066903.98 4081.25 | 46.40
25. -2066748.43 4170.60 | -2066526.69 4001.82 | 59.12
26. -2066067.79 2292.29 | -2065860.31 2094.70 | 46.35
27. -2065336.69 2107.87 | -2065125.92 1900.40 | 43.07
28. -2064605.58 1913.57 | -2064378.35 1692.92 | 26.16
29. -2064199.15 3558.56 | -2063961.13 3356.40 | 47.43
30. -2037464.39 6512.45 | -2037075.56 5864.40 | 50.66
31. -2036721.82 6825.68 | -2036338.39 6199.23 | 22.54
32. -2036875.74 6366.64 | -2036497.71 5742.89 | 21.95
33. -2040224.67 7150.18 | -2039706.23 6575.03 | 161.54
34. -2041064.45 7144.78 | -2040732.32 6569.63 | 26.74
35. -2044701.68 -15548.03 | -2044564.34 -16024.90 | 23.64
36. -2043992.10 -15723.52 | -2043824.24 -16223.28 | 47.60
37. -2043789.90 -14907.12 | -2043610.60 -15383.99 | 34.35
38. -2040615.94 -14820.44 | -2040453.30 -15349.97 | 20.77
39. -2039594.70 -15081.43 | -2039485.02 -15588.26 | 47.85
Number of points: 39
Residual mean average: 57.082951
... Continuing from the initial example session in R ...
An affine transformation is based on a linear mapping between two coordinate-spaces. Testing for non-linearity (i.e. higher order model terms) can be a useful diagnostic in choosing the simpler affine transform.
Compute the difference between the good and bad coordinates
Generate two linear models:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 302.38 3.97 76.174 < 2e-16 ***
poly(nx, 3)1 1091.17 35.82 30.466 < 2e-16 ***
poly(nx, 3)2 165.59 32.78 5.051 1.71e-05 ***
poly(nx, 3)3 -51.21 28.74 -1.782 0.0843 .
poly(ny, 3)1 417.35 31.97 13.056 2.30e-14 ***
poly(ny, 3)2 -18.50 30.04 -0.616 0.5425
poly(ny, 3)3 19.41 35.49 0.547 0.5882
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Is one model significantly better than the other?
Analysis of Variance Table Model 1: sqdiff ~ nx + ny Model 2: sqdiff ~ poly(nx, 3) + poly(ny, 3) Res.Df RSS Df Sum of Sq F Pr(>F) 1 36 52093 2 32 19665 4 32428 13.192 1.865e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Check visually:
Testing for linearity: Two visualizations of the deviance between coordinates positions, at the control point locations.
Conclusions:
It seems that a second order term was only warranted along the x-direction. The more complex model based on 3rd-order polynomials results in a significantly lower RMSE (about 10 meters lower), and is shown to be a better descriptor of variance in the test of nested models.
At the map scale in which the corrected data will be presented, the extra accuracy suggested by the more complex model (coordinate transformation function) is negligible. This allows for the simpler model, which can be directly used by the convenient ST_Affine() function in PostGIS for the heavy-lifting.
The following is a proof of concept query showing how a PostgreSQL query could give the transformation parameters required for an affine transformation. The source of the procedure is from GRASS at transform.c.
Note that the v.transform code expects the input file (control points) to be in the form ax ay bx by, where 'a' is the starting coordinate system (the bad coordinates in the previous example) and 'b' is the target system (the good coordinates). The
SQL version is as follows:
Sample Session
This query requires a table called link with the following fields - gid (primary key), a_x, a_y , b_x, b_y. The 'a' values are the 'from' coordinates and the 'b' values are the 'to' coordinates. Using the attached control points the result of this query is:
b0 | b1 | b2 | b3 | b4 | b5
------------------+---------------------+-------------------+------------------+------------------+---------------------
-28138.394850347 | -0.0135202854235867 | 0.997400773420259 | 5017.08164289594 | 1.00231638907948 | 0.00918961946271679
These are the exact results from 'R'
Based on the results of the proof of concept example developed previously, a single function was developed in the procedural language for the PostgreSQL database system called PL/pgSQL. The only input parameter to this procedure is a text string that results in a table in the following format:
Note the table must have the above fields but they can be in any order and can have additional fields. The gid field must be unique for each record.
The use of a SQL query found to be the simplest way to avoid the difficulties in programing the procedure without the need for temporary tables. An added benefit is the control point data can be in almost any format as long as it can be arranged in the format specified above. For example if the "from" points are in a geometry column (the_geom) in a table called from_pts and the corresponding "to" points are in a similar table called to_pts with a common attribute called "link_id", then the query would be:
Other table layouts and queries are possible depending on the manner in which the control points are collected.
The following is the SQL code to add a new procedure called trans_param() into a PostGIS database:
Currently there is no error checking in the code if the determinant is zero.
To use the procedure simply use: SELECT trans_param('my SQL text')
The data from this series of articles is stored in a table called link. The id of the points is gid and the "from" values are b_x and b_y. The "to" values are a_x and a_y. Thus the query is:
The result of the procedure is:
0.997546509279282;-0.00917177514909895;0.0135300872142122;1.00246938174737;-5301.39933295548;28155.8822879205
Which matches the results returned from GRASS and R. Additional queries will be developed to give a table of residuals, the RMS error and the actual transformation of geometry.
This is a start for discussions to create a series of function to perform an affine transformation of a PostGIS data set using a table of control points.
The first step is what is the format or layout of the control points?
I don't think point geometries is necessarily a good idea. The points must have an exact 1 to 1 relationship. To assure this you either have to maintain absolute integrity on the keys between two tables or have two geometries in one table. Both would be a hassle. An interleaved table format also would give me problems because the queries get really difficult and what if one version used "good" points in the odd rows and the next one put them in the even ones? In addition, if you use a seguential id and add one row "good" then add a "bad" one, then realize the bad one is really bad and delete it then the gid's will not be odd/even anymore in relation to "good" vs "bad". I think ESRI has the right idea here (God did I actually say that?). Look at the link table interface when you are georeferencing an image. It gives the id, XSource, YScource, Xmap, YMap, and the residual error. This could easily be done with the code - just one more query.
Given a table in this format a query could be done to give the RMS error. Also the table could be returned with the residual errors. When the user is happy then the table and the geometry could be input to transform the geometry.
Bruce
News
It looks like the NCSS has constructed a web-based, SQL interface to their main database. This new tool was recently highlighted in issue 40 of the NCSS Newsletter and looks like a promising new tool with good documentation. A listing "web services" offered by the NCSS is listed on this page.
Overview
The analysis of SSURGO data is complicated by a series of hierarchical, one-to-many tabular relationships between spatial and attribute data. Understanding the structure of the SSURGO database is critical for correct interpretation and aggregation of soil properties. Before undertaking any SSURGO-based analysis please take some time to become familliar with the details on the product, including intended uses, minimum mapping units sizes, and other important details. In addition, becoming familliar with the SSURGO metadata is a must. The metadata page includes detailed descriptions of table structure, column names and units, as well as important information on the sources of much of the tabular data included in a SSURGO database.
SSURGO data can be downloaded by survey area from The Soil Datamart, with spatial data delivered in shapefile format and attribute data delivered as plain text. Unfortunately, an M.S. Access database template is required to utilize SSURGO attribute data as delivered from the Soil Datamart. For assistance with this procedure please see the NRCS SSURGO page. Using this approach, most analysis of SSURGO must be done survey-by-survey within a GIS environment. For general instructions see this document.
Several online applications allow for simple interaction with the SSURGO database without the need for a GIS or RDBMS knowledge. Some examples include:
An Open Source Approach to SSURGO
We have developed an alternate approach to working with SSURGO data using PostGIS, a spatially-enabled version of the popular and open source relational database management system Postgresql, to store all spatial and tabular data for 138 survey areas. This approach facillitates rapid access, analysis, and aggregation of over half a million soil polygons. SQL (structured query language) is used to directly interact with both soil spatial and attribute data. If other forms of spatial data are imported into PostGIS (such as landcover, climatic data, etc.) nearly all spatial and attribute analysis can be done entirely from PostGIS. A series of examples illustrating common tasks will be presented in the following pages.
General Approach to Working with SSURGO (also outlined in this document) (more ideas on map unit composition)
SSURGO discontinuity example: Boundary between Glenn and Colusa counties, illustrating differences in soil survey vintage.
Notes on the Format of SSURGO
SSURGO Database Structure Diagram
Building and Saving SQL code snippets
A simple text editor is the best environment for working on (and saving!) your SQL queries. Linux users are encouraged to use either 'Kate' or 'Kwrite'. Windows users should look into Notetab Light or SciTE. Mac users should check out TextWrangler.
Connecting to the database
In general, the simplest way to interact with our composite soil survey database is by connecting to the best with an SSH client. Mac/Linux users have this functionality built-in. Windows users will need to use something like Putty, or Xming. Once an ssh connection with the beast has been setup, you can connect to the database with the following:
where psql is the postgresql client program, -U soil means connect as the user called "soil", and ssurgo_combined is the name of the actual database.
You can quit from the database shell using \q (followed by enter). Typing \? (followed by enter) will give a list of commands available in the psql shell.
Query structure and SQL
Numerous resources exisit for learning about SQL. See the attached PDF presentation at the bottom of this page for some SSURGO-related examples. I would recommend looking at the first couple chapters from the PostgreSQL book by Douglas and Douglas (the big purple book on the shelf). For an interactive learning approach SQL Zoo seems like a good start. In general most queries will have the format:
A simple query: column selection and filtering
Once connected it is possible to issue queries to the database. The results of a query are normally returned to the screen. For example, asking for the horizon boundaries and horizon names from the componenent identified with the given cokey (467038:646635) might look like this:
A more complicated query: column calculation and aggregation
Compute the total water holding capacity for a given component, identified by (467038:646635). Note that several operations are being performed on the data:
-- the query select cokey, sum( (hzdepb_r - hzdept_r) * awc_r) as component_whc, sum( (hzdepb_r - hzdept_r) ) as depth from chorizon where cokey = '467038:646635' group by cokey ; -- the results cokey | component_whc | depth ---------------+---------------+------- 467038:646635 | 27.91 | 183
Note that it is also possible to save the results of a query into a new table. This is the usual strategy when performing any analysis that returns geometry (soil polygons, etc.) as geometry data cannot be visualized in a text display. Tables created in this manner can be exported from the database for map creation or further analysis. Spatial tables can be viewed directly by applications like QGIS or Mapserver. A simple method of accessing PostGIS data is not yet possible with ArcGIS.
Just Checking
-- NAD27 to NAD83 echo 119d7\'4\"W 36d23\'13\"N | cs2cs +proj=latlong +datum=NAD27 +to +proj=latlong +datum=NAD83 -f "%.6f" -119.118718 36.386894 0.000037
Selection, Filtering and Sorting
The resulting table can be copied to a CSV file (in the current working directory) like this: \copy yolo_comp TO 'yolo_comp.csv' CSV HEADER
If you are done with the table, remove it with the following SQL: DROP TABLE yolo_comp ;
Overview:
Two methods for the selection of the largest components (based on the comppct_r column) within map units. This approach to selecting a single component per map unit is appropriate in situations where a single representative feature or property is sought. The partitioning of components within a map unit is closely related to the map unit type: complex, association, consociation or an "undifferentiated group". A breakdown of the number of components per each map unit type is summarized by the following query (26400 map units / 45971 major components):
Method 1: Filtering component percentages with the DISTINCT keyword.
Finding the largest component (based on comppct_r) per mukey
SELECT count(mukey) FROM mapunit WHERE areasymbol = 'ca113';
Overview
This example illustrates a component-percent weighted mean water storage, based on total profile storage values. The comppct_r column is used as a weight to adjust the influence of the profile water storage values for each component: larger weights (larger components) result in more influence. When performing weighted means, use care in the selection of an appropriate weighting parameter. A selection of a weighting variable not direclty related to the property which is being averaged can result in an effect known as Simpson's paradox.
Some background on Soil Water Holding Capacity and Irrigation Management.
The general formula for computing profile water storage is defined as follows. Note that the awc_r column in the SSURGO database (chorizon table) is pre-adjusted to compensate for coarse fragments (see relevant SSURGO metadata below).
Calculate the water storage within a given component
component_whc = profile_sum_awc = sum( (hzdepb_r - hzdept_r) * awc_r )
Calculate the weighted mean value of the profile water storage, for a given map unit
mapunit_whc = sum( (comppct_r/100) * component_whc ) / sum( comppct_r/100 )
From the SSURGO 2.1 table definitions:
Create indexes and register the new geometry:
Yolo County Water Holding Capacity Map: Profile total water holding capacity as computed by component percentage weighted average.
ogr2ogr ywhc.shp PG:"dbname=ssurgo_combined user=xxxx password=xxxx" yolo_whc
Premise
Compute a series of weighted-average soil texture fractions (sand, silt, clay), for every component, of every map unit in Yolo County. These values will be further weighted by the spatial distribution of each map unit.
Simple Visualization with R
Yolo County Common Soil Textures: Blue symbols mark common soil textures and their relative size denotes abundance. The orange symbol marks the weighted average soil texture for all of Yolo County.
Overview
A simple association between dated landforms and soil series name [1] was used to extract soil polygons from a composite soil survey database.
| Soil Series | Associated Formation | Approximate Age (1000 yrs ago) |
| Redding | Laguna | 1600 - 2000 kya |
| Corning | Laguna | 1600 - 2000 kya |
| Keyes | Laguna | 1600 - 2000 kya |
| Whitney | Turlock Lake | 500 - 700 kya |
| Montpellier | Turlock Lake | 500 - 700 kya |
| Rocklin | Turlock Lake | 500 - 700 kya |
| Snelling | Riverbank | 100 - 300 kya |
| San Joaquin | Riverbank | 100 - 300 kya |
| Exiter | Riverbank | 100 - 300 kya |
| Madera | Riverbank | 100 - 300 kya |
| Hanford | Modesto | 10 - 40 kya |
| Grangeville | Holocene | < 10 kya |
| area_ac | formation |
|---|---|
| 295868 | Riverbank |
| 253424 | Modesto |
| 151149 | Turlock Lake |
| 121085 | Laguna |
| 96981 | Holocene |
(5 rows)
ogr2ogr east_side_all.shp PG:"dbname=ssurgo_combined user=xxxx password=xxxx host=xxxx" east_side_all
References:
Example 1: The location of seasonaly wet soils via two methods: hydricrating and USDA Soil Taxonomy interpretation.
Example 1.1 Extract a list of wet components, and sum area based on component (series) name
musym | mukey | muname | wet_area_ac -------+--------+----------------------------------------------------------------+------------- Sc | 459268 | Sacramento clay | 35710.90 Mf | 459244 | Marvin silty clay loam | 18877.18 Sg | 459272 | Sacramento soils, flooded | 12273.51 Cn | 459219 | Clear Lake soils, flooded | 11665.84 Cc | 459216 | Capay soils, flooded | 11029.83 Sv | 459288 | Sycamore complex, drained | 8410.41 St | 459286 | Sycamore silty clay loam, drained | 6900.77 Ck | 459218 | Clear Lake clay | 6737.88 Sa | 459266 | Sacramento silty clay loam | 6014.48 Sp | 459283 | Sycamore silt loam, drained | 5537.18 Sw | 459289 | Sycamore complex, flooded | 5041.96 Wb | 459301 | Willows clay | 4950.34 Ss | 459285 | Sycamore silty clay loam | 4859.59 Pb | 459254 | Pescadero silty clay, saline-alkali | 4700.77 So | 459282 | Sycamore silt loam | 3953.42 Rh | 459262 | Riverwash | 3698.88 Pc | 459255 | Pescadero soils, flooded | 3589.01 Tb | 459292 | Tyndall very fine sandy loam | 3299.03 Ob | 459252 | Omni silty clay | 3174.44 [...]
Example 2: Classify soils according to shrink-swell capacity of the top 1 meter of soil, weighted by horizon thickness and component percent.
musym | muname | muacres | lep | lep_class -------+------------------------------------------------------------------+---------+------+----------- Ya | Yolo silt loam | 39698 | 2.52 | Low Sc | Sacramento clay | 34886 | 7.50 | High Ca | Capay silty clay | 33465 | 7.50 | High MrG2 | Millsholm rocky loam, 15 to 75 percent slopes, eroded | 30118 | 1.50 | Low Rg | Rincon silty clay loam | 24580 | 6.36 | High BrA | Brentwood silty clay loam, 0 to 2 percent slopes | 23045 | 7.50 | High CtD2 | Corning gravelly loam, 2 to 15 percent slopes, eroded | 22080 | 5.34 | Moderate Mf | Marvin silty clay loam | 20970 | 6.60 | High DaF2 | Dibble clay loam, 30 to 50 percent slopes, eroded | 18612 | 7.11 | High SmE2 | Sehorn-Balcom complex, 15 to 30 percent slopes, eroded | 17794 | 6.17 | High TaA | Tehama loam, 0 to 2 percent slopes | 16622 | 3.75 | Moderate BdF2 | Balcom-Dibble complex, 30 to 50 percent slopes, eroded | 16405 | 5.73 | Moderate SmD | Sehorn-Balcom complex, 2 to 15 percent slopes | 16117 | 6.50 | High BaF2 | Balcom silty clay loam, 30 to 50 percent slopes, eroded | 12637 | 4.50 | Moderate Sg | Sacramento soils, flooded | 12258 | 6.27 | High Cn | Clear Lake soils, flooded | 11666 | 6.92 | High SmF2 | Sehorn-Balcom complex, 30 to 50 percent slopes, eroded | 11226 | 6.33 | High Cc | Capay soils, flooded | 11030 | 7.50 | High Sv | Sycamore complex, drained | 9241 | 4.18 | Moderate Ms | Myers clay | 8938 | 7.50 | High PfF2 | Positas gravelly loam, 30 to 50 percent slopes, eroded | 7920 | 5.34 | Moderate St | Sycamore silty clay loam, drained | 7839 | 4.50 | Moderate Ck | Clear Lake clay | 6946 | 7.50 | High Ra | Reiff very fine sandy loam | 6847 | 1.50 | Low [...]
http://www.rbt.ca/autodoc/index.html postgresql_autodoc -d db_name see dot output. dot -Tpng -o output.png db_name.dot
Data points collected in the field via GPS, or in the office via maps can often come from a variety of projections and or datum/ellipsoids. The PROJ tools and libraries can be used to perform forward and reverse projection on lists of points stored in a text file. More information on PROJ can be found here: http://trac.osgeo.org/proj/ . Many applications that perform forward and inverse projection operations rely on the Proj library for the actual transformatoin of coordinates. Some examples include GDAL/OGR and GRASS, which provide a convenient interface for converting large datasets from one coordinate system to another. For problems see the FAQ. Tests confirm that datum transforms within the continental US are nearly identical with results from NADCON.
Command names are in boldface type, input files are red, and output files are in green.
Forward and Inverse Projections
Conversions Between Coordinate Systems
Overview
Operations on geographic data are most efficient when the input files have identical spatial parameters: i.e. coordinates system and datum / ellispoid. The GDAL/OGR tools and libraries can be used to "reproject" downloaded geographic data from one coordinate system to another, both for vector and raster data. More information on this software can be found here: http://www.remotesensing.org/gdal/index.html The following examples illustrate some of the ways in which the GDAL/OGR tools can be used to manipulate the spatial parameters of some common types of geographic data: SSURGO (ESRI Shapefile format) and aerial imagery (GeoTiff format).
How Does it Compare with ArcToolbox?
A SSURGO boundary polygon was used in a quick test to see how similar the results of a vector re-projection done by ArcMap and OGR were. The input vector (GCS NAD83) was projected to UTM zone 11 NAD83 by both OGR and ArcMap 9.0, and then imported into GRASS. Both vectors were then exported as ascii, for simple access to the line vertex coordinates. A simple AWK script was used to numerically compare the differences in both X and Y coordinates:
paste arcmap.ascii ogr.ascii | \
awk '
/\./{print ($1-$3), ($2-$4) ; x+=($1-$3) ; y+=($2-$4); n++ }
END{print "avg dx: "(x/n), "avg dy: "(y/n), "n: "n}
'
The resulting average dx and dy values were: dx: 2.98838e-06 m and dy: -1.23948e-07 m with 5035 samples: looks good to me. Next time, a datum shift may be more telling of any subtle differences to be aware of. Update: NAD83 -> NAD27 Datum shift experiment looks like datum shifts are comparible as well.
Examples:
Command names are in boldface type, input files are red, and output files are in green. Figures are provided for additional context, and are referenced by input/output file name in the example code.
Vector Operations
Raster Operations [raster tutorial]
Reading MrSid format files details
Comment out this line, and re-compile. Thanks to Mateusz Loskot on the GDAL IRC channel for pointing this out.
About R
R is a general-purpose, command-line based, environment for working with data. R is based on a functional approach to working with vectors and matrices of data. R is a convenient environment for processing, analyzing, and plotting data. The statistical tests and model building facilities of R are rivaled by no other single application.
R in the News
Getting Started
Searching for Information
R with Geographic Data
Misc. Articles
Overview
Perform some temporal aggregation (by day and by week) of the amount of data entry completed in Postgresql, and plot the results in R. See resulting figure at the bottom of the page. Note that this requires the Rdbi and RdbiPgSQL packages. Hints on installing these packages can be found on this page...
Weekly Aggregation hints from the psql manual page
Daily Aggregation hints from the psql manual page
R Example

Decagon Sensors: EC-5 (moisture) and ECT (temperature)
Premise
Simple demonstration of working with time-series data collected from Decagon Devices soil moisture and temperature sensors. These sensors were installed in a potted plant, that was semi-regularly watered, and data were collected for about 80 days on an hourly basis. Several basic operations in R are demonstrated:
Process the raw sensor values with standard calibrations
Decompose each time series into additive components
Additive Time Series Decomposition: Temperature
Additive Time Series Decomposition: Volumetric Water Content
Auto-Correlation Function (ACF)
Soil Moisture and Temperature ACF: Auto-correlation function of each time series, and cross-ACF.
Interesting Results
Variation in temperature with time dominated by diurnal fluctuations superposed over underlying fluctuations caused by building heating/cooling system. The magnitude of the diurnal cycle appears to be related to the moisture content- as expected due to high heat capacity of water. Diurnal variation in moisture values appears to account for less than < 2% absolute change in volumetric water content.
Examples based on a random data set (see example code below), illustrating some of the differences between the K-means and C-means clustering methods as implemented in R. Next time an example with soil profile data collected from the Pinnacles National Monument soil survey efforts. An online version of the PINN soil survey will be available soon here.
Articles:
Example in R:
Sample functions and ideas for accessing the R built-in colors. Further examples on converting soil colors to RGB triplets, or for the selection of optimal colors for a thematic map please see the examples linked at the bottom of this page. An excellent discussion on the use of color for presenting scientific data is presented in this paper by Zeileis, Achim and Hornik, Kurt.
R Color Selection: Simple figure illustrating the layout() function to create a plot of the built-in R colors palettes.
Simple Color Display
A Queryable color picker (as suggested by Gabor Grothendieck on the R-help mailing list)
Setup the plot layout, and plot both examples
The Munsell color system was designed as a series of discrete color chips which closely approximation to the color sensitivity of the human eye. The description of color via three variables tied to perceptible properties (hue, value, and chroma) under a standardized illuminant (sunlight on a clear day) makes the Munsell system a good choice for recording and interpreting soil color data. However, numerical analysis of colors encoded in the Munsell system is difficult because they are from a discrete set of color chips and referenced by values that include both letters and numbers. Rossel et. al. (2006) give a good background on various color space models and their relative usefulness in the realm of soil science. The conversion of Munsell soil colors to RGB triplets, suitable for displaying on a computer screen or printing, is made complicated by the numerous operations involved in converting between color spaces. Figure 1 shows all possible (both real and unreal) Munsell color chips in the L*U*V color space. Figure 2 shows some of the common soil color chips in the same color space. Figures 2 through 5 depict common soil colors in the RGB color space, visualized both in R and POVRAY. Example R code on the conversion is given below.
Munsell color data can be downloaded here.
Color conversion equations here.
References:
Setup environment and load lookup table data
Convert xyY to XYZ [Equation Reference]
Perform Chromatic Adaption Functions in the colorspace package, and sRGB profiles assume a D65 illuminant [Reference]
Convert XYZ (D65) to sRGB (D65), step 1 this assumes that XYZ is scaled to (0,1) [Reference Primaries for sRGB]
Convert XYZ (D65) to sRGB (D65), step 2 (sRGB, gamma = 2.4) [Conversion Function to sRGB]
Choosing the right colors for classes in a thematic map can be a difficult task. The ColorBrewer website provides an interactive tool for browsing numerous color combinations. Each of the color combinations presented on the ColorBrewer website are the culmination of numerous color interpretation studies. In addition, there is a list of special color combinations suitible for audiences which may include color blind individuals.
The R package RColorBrewer adds the color brewer color combinations as well as functions for generating new combinations to R. Figure 1 demonstrates the available color combinations, as returned by the function display.brewer.all.
An example R session:
When the relationship between two variable is (potentially) dependent on a third, categorical variable ANCOVA (analysis of covariance), or some variant, is commonly used. There are several approaches to testing for differences in slope/intercepts (in the case of a simple linear model) between levels of the stratifying variable. In R the following formula notation is usually used to test for interaction between levels of a factor (f) and the relationship between two continuous variables x and y: y ~ x * f. A simple graphical exploration of this type of model can be done through examination of confidence intervals computed for slope and intercept terms, for each level of our grouping factor (f). An example of a fictitious dataset is presented below. Note that this a rough approximation for testing differences in slope/intercept within a multi-level model. A more robust approach would take into account that we are trying to make several pair-wise comparisons, i.e. something akin to Tukey's HSD. Something like this can be done with the multcomp package. For any real data set you should always consult a real statistician.
Example Multi-Level Model: each panel represents a model fit to y ~ x, for group f
Example Multi-Level Data
Example Multi-Level Model: Confidence Intervals: parameter estimates along with 95% confidence interval, within each level of our grouping factor (f).
Automated Plotting of Parameter Confidence Intervals
The first two lines in the output below are testing the hypothesis that the slope and intercept term for level 'a' are not different than 0. Subsequent hypothesis tests are relative to the first 'level' in the dataset. In this case we are testing the hypothesis that intercept and slope terms for levels 'b' through 'j' are not different than the corresponding terms for level 'a'. From the output below we can see that none of the intercept terms (levels 'b' through 'j') are different than for 'a', and that the slope term for level 'd' is only marginally "different" (p=0.0625) than the slope term for 'a'.
Testing Model Terms
Small update to a similar thread from last week, on the comparison of slope and intercept terms fit to a multi-level model. I finally figured out (thanks R-Help mailing list!) how to efficiently use contrasts in R. The C() function can be called within a model formula, to reset the base level of an un-ordered factor. The UCLA Stats Library has an extensive description of this topic here. This approach can be used to sequentially test for differences between slope and intercept terms from a multi-level model, by re-setting the base level of a factor. See example data and figure below.
Note that the multcomp package has a much more robust approach to this type of operation. Details below.
Example Multi-Level Data
Default Contrasts (contr.treatment for regular factors, contr.poly for ordered factors)
Setting the "base level" in the Model Formula This allows us to compare all slope and intercept terms to the slope and intercept from level 4 of our factor ('d' in our example).
Testing with Multcomp Package using data from above example
Panel function for visualizing univariate statistics
Panel function for printing joint correlation statistic
Example usage with built-in datasets
Overview
Sample application of the RGL package. This package allows for the creation of interactive, 3D figures, complete with lighting and material effects. Try demo(rgl) for an idea of what is possible.
A random number generator sphere (RNG sphere) was created based on the suggestions in Keys to Infinity by Clifford A. Pickover, pp. 237-239. The RNG sphere can be used to test the robustness of a random number generator. Three random number generators were tested: runif() from R, rand from Excel, and a logistic-derived psudo-random number generator. The location (x,y,z) and color of the spheres are based on the sequence of random numbers (Pickover, 1995). An ideal RNG shpere should have no discernable patterns. Note that the logistic-derived random numbers show distinct correlation in the RNG sphere. Excel random number list, and source code (R) are attached at the botom of the page.
RGL sample application: 3d interactive interface to a random number generator sphere. Random numbers from runif() function in R.
RGL sample application 2: Excel random number visualization: 3d interactive interface to a random number generator sphere. Random numbers from rand() function in MS Excel.
RGL sample application 2: Random numbers from the logistic function: 3d interactive interface to a random number generator sphere. Random numbers from the logistic function (see notes), implemented in R.
Random Number Generator (RNG) Sphere Function Definition
Sample
Generate some data
The dataframe 'soil' looks like this:
depth prop 1 0 13.80257 ** note that these are the same 2 12 13.80257 ** note that these are the same 3 24 18.40298 4 36 13.37446 5 48 13.27973 6 60 14.65288 7 72 16.07339 8 84 15.97451 9 96 16.29970 10 108 16.32155 11 120 14.63699 12 132 13.26486 13 144 13.81730
Plot the data:
Notes:
Helper Function (copy this into your R session first)
Load Data and Packages
Example 1
Depth Profile Example 1: sand, silt, and clay vs. depth for three pedons
Example 2
Premise:
Some examples on how to prepare and present data collected from an XRD analysis. The clay fraction from seven horizons was analyzed by XRD, using the five common treatments: potassium saturation (K), potassium saturation heated to 350 Deg C (K 350), potassium saturation heated to 550 Deg C (K 550), magnesium saturation (Mg), and magnesium + glycerin saturation (Mg+GLY). These data files have been attached, and can be found near the bottom of the page.
Plotting the entire data set with lattice graphics:
Example XRD plot with lattice graphics: 7 horizons and 5 treatments
Locating relevant peaks in an X-ray diffractogram is an important step in identifying phyllosilicate minerals in soils. An automated approach to finding peaks in any dataset was presented by Martin Maechler, contributed to the R-Help mailing list Nov 25, 2005. paste these functions into an R session to use them
Commands to find and plot the peaks, based on suggestions by peaks() function author.