Well, I made some mistakes in the last post, not the least of which is I used the wrong flag for creating an attribute field with elevation. What follows is a little more sophisticated. It takes us from a series of DEM tiles from which I want 2-foot and 5-foot contours (using gdal_contour), and then dumps those shapefiles into PostgreSQL using shp2pgsql.
First we prep the database. Again, I used pre-built binaries from Kyng Chaos.
I created a database called “Vinton” using my favorite PostgreSQL management tool, PGAdminIII (actually, it’s the only one I’ve ever tried…).
So, let’s make Vinton spatially aware:
sudo su - postgres -c '/usr/local/pgsql/bin/createlang plpgsql Vinton'
sudo su - postgres -c '/usr/local/pgsql/bin/psql -d Vinton -f /usr/local/pgsql/share/contrib/postgis.sql'
sudo su - postgres -c '/usr/local/pgsql/bin/psql -d Vinton -f /usr/local/pgsql/share/contrib/spatial_ref_sys.sql'
Vinton county is in Ohio State Plane South. Our dataset is a foot version of the State Plane South HARN. Looking it up on SpatialReference.org, I get a description not unlike the PostGIS insert statement below, but I changed the EPSG number from 102323 to 1023236 (arbitrarily), and changed the units to feet (from the default for HARN, which is meters):
UNIT["Foot_US",0.30480060960121924]
so that the full insert is:
INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 91023236, 'esri', 1023236, '+proj=lcc +lat_1=38.73333333333333 +lat_2=40.03333333333333 +lat_0=38 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +units=m +no_defs ', 'PROJCS["NAD_1983_HARN_StatePlane_Ohio_South_FIPS_3402",GEOGCS["GCS_North_American_1983_HARN",DATUM["NAD83_High_Accuracy_Regional_Network",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",600000],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-82.5],PARAMETER["Standard_Parallel_1",38.73333333333333],PARAMETER["Standard_Parallel_2",40.03333333333333],PARAMETER["Latitude_Of_Origin",38],UNIT["Foot_US",0.30480060960121924],AUTHORITY["EPSG","1023236"]]');
Now let’s create some tables to house our future contour data:
CREATE TABLE contours_2
(
gid serial NOT NULL,
id integer,
elev numeric,
the_geom geometry,
CONSTRAINT contours_2_pkey PRIMARY KEY (gid),
CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL),
CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = 91023236)
)
WITH (OIDS=FALSE);
ALTER TABLE contours_2 OWNER TO postgres;
CREATE TABLE contours_5
(
gid serial NOT NULL,
id integer,
elev numeric,
the_geom geometry,
CONSTRAINT contours_5_pkey PRIMARY KEY (gid),
CONSTRAINT enforce_dims_the_geom CHECK (st_ndims(the_geom) = 2),
CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = 'MULTILINESTRING'::text OR the_geom IS NULL),
CONSTRAINT enforce_srid_the_geom CHECK (st_srid(the_geom) = 91023236)
)
WITH (OIDS=FALSE);
ALTER TABLE contours_5 OWNER TO postgres;
And finally, we can loop through all of our files, and dump the results first to a shapefile, then pipe that from shp2pgsql directly into the database:
#!/bin/bash
x=0
for f in $( ls *.txt); do
x=`expr $x + 1`
echo $x $f
echo Two Foot:
gdal_contour -i 2 -a elev $f $f.shp
/usr/local/pgsql/bin/shp2pgsql -s 91023236 -a $f.shp contours_2 | /usr/local/pgsql/bin/psql -U postgres -d Vinton
echo Five Foot:
gdal_contour -i 5 -a elev $f 5foot/$f.shp
/usr/local/pgsql/bin/shp2pgsql -s 91023236 -a 5foot/$f.shp contours_5 | /usr/local/pgsql/bin/psql -U postgres -d Vinton
done
Now, I’ll probably convert all those numeric contour fields to integer with a cast statement, but tonight, I let the little laptop do the work… .
