Simple Python Interface to USGS TNM Elevation Service

 
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.

 
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

Text processing with the UNIX toolbox

 
Working with DBF files

• Extract the text contents of a DBF file with dbfxtrct, select out the 3rd and 2nd columns with awk, keep only unique rows with sort and uniq, and finally saving to a new file called 'ca_fips-name_list.txt' with shell redirection.
  
  # dbfxtrct utility
  dbfxtrct -ica_counties_aea83.dbf | awk -F"," '{print $3","$2}' | sort -g | uniq > ca_fips-name_list.txt
• Other DBF tools include:
  • dbf2mysql
  • dbf2pg
  • dbfadd
  • dbfcreate
  • dbfdump
  • dbfutil1
  • dbfxtrct

Customizing BibTeX

 
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:

• when citing named chapters in a book, use the incollection entry type
• when using makebst get the style as close as you can to the target, then tweak the resulting .bst file

 
Documentation:

http://www.andy-roberts.net/misc/latex/latextutorial3.html

 
A Custom Bib Style:

Latex to ODT conversion tips

 
Conversion of a .tex document to open document format (RTF)

 
Conversion of a .tex document to open document format (ODT)

 
Some notes:

1. disable hyperref package to prevent odd bug in citation listings
2. input document should be latex-safe, i.e. eps figures
3. need a working java installation
4. commands like \onehalfspacing or begin{singlespace}...end{singlespace} may not work (setspace package)

 
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.

1. Some related pages
   • http://www.physik.tu-cottbus.de/~george/unix/emf4unix.html
   • http://boshoff.za.net/linux/archives/2003_08.html
2. download, hack, and compile the libEMF code (details)
3. download, configure, and install pstoedit, note that you may have to manually define where the libEMF header and lib files are:
   
   ./configure --with-libemf-include=/usr/local/include/libEMF/ \
   --with-libemf-lib=/usr/local/lib/
4. run ldconfig as root to refresh the library cache
5. [Finally] Convert a PDF file into an OpenOffice-compatible EMF file (well sort of...)
   
   pstoedit -pta -f emf:"-OO -p -drawbb" file.pdf file.emf
6. It seems like figures are clipped to smaller-than-expected bounding boxes unless the special flag -drawbb is used. One side-effect of this however is a box drawn around your figure. Also, the -pta flag is needed to ensure quality text placement. Even though the resulting EMF files may look odd on a Linux machine, they appear to work as expected in Windows.



Making better looking tables with \ctable

 
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}

Importing and Exporting

Tips

• Attribute data stored with a shapefile (the .dbf file) can only use column names shorter than 14 characters. If you try and export a postGIS table that has column names longer than 14 characters, they will be truncated when converted to a shapefile.
• If a column in postGIS is defined as numeric, ogr2ogr will not correctly interpret the values in this column: i.e. floating point values less than 1 will be truncated to 0. In order to avoid this problem always define your colums as double precision or float. For dynamically generated tables, cast affected columns to a defined numeric precision: select colum_a::numeric(7,3) where '7' refers to the total number of digits in the column, and '3' refers to the number of decimal places.

Geometry and attribute data to GIS file format

 
GDAL/OGR tools This approach allows simultaneous conversion of coordinate systems, but is less flexible with respect to generation of new tables in PostGIS.

• Import
   
  ogr2ogr -f "PostgreSQL" PG:'dbname=ssurgo_combined user=xxxx password=xxxx host=postgis.server.edu' input_file.shp
• Export
   
  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.

• Import
   
  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.

• Export
   
  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.

Attribute data to text format

 
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

(10 rows)

Example Spatial SQL Operations on Point, Line and Polygon Geometry Types

• Points
  • Example
  • Example
  • Example
• Lines
  • Example
  • Example
  • Example
• Polygons
  • Example
  • Example
  • Example
• Mixed
  • Example
  • Example
  • Example

Affine Transformation Operations in PostGIS

 
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).

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

          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

Transformation Matrix
| xoff a b |
| yoff d e |
-------------------------------------------
5301.399323 1.002469 0.009172 
-28155.882288 -0.013530 0.997547 
-------------------------------------------

               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

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 

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

        b0        |         b1          |        b2         |        b3        |        b4        |         b5          
------------------+---------------------+-------------------+------------------+------------------+---------------------
 -28138.394850347 | -0.0135202854235867 | 0.997400773420259 | 5017.08164289594 | 1.00231638907948 | 0.00918961946271679

0.997546509279282;-0.00917177514909895;0.0135300872142122;1.00246938174737;-5301.39933295548;28155.8822879205

Analysis of SSURGO Data in PostGIS: An Overview

 
Note
Several of the pages documenting soil survey related queries have been taken down for a couple of days while I update the syntax.

 
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 assist