Import script for NHDPlus.

nhdplus_import.sh

#!/bin/sh
#
# Create a database for NHDPlus to live in and load it.
#
# This presumes that the active user is a (Postgres) super-user and that all
# dependencies are installed.

# set -e -x
set -e

DBNAME=nhdplus
SRID=4269
ENCODING=LATIN1
FLOWLINE_TABLE=flowline
AREA_TABLE=area
WB_TABLE=waterbody
PT_TABLE=point
LINE_TABLE=line

createdb $DBNAME
psql -d $DBNAME -c "CREATE EXTENSION postgis;"

# for f in *.7z; do
#   7z x $f
# done

## Flowlines
flowlines=$(find . -iname nhdflowline.shp)

set -- $flowlines
echo "Creating $FLOWLINE_TABLE schema"
shp2pgsql -p -D -t 2D -s $SRID $1 $FLOWLINE_TABLE | psql -d $DBNAME -q

for f in $flowlines; do
  echo "Importing $f..."
  shp2pgsql -a -D -t 2D -s $SRID -W $ENCODING $f $FLOWLINE_TABLE | pv | psql -d $DBNAME -q
done

psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${FLOWLINE_TABLE}_geom_gist ON $FLOWLINE_TABLE USING GIST(geom);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${FLOWLINE_TABLE}_ftype_idx ON $FLOWLINE_TABLE(ftype);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${FLOWLINE_TABLE}_comid_idx ON $FLOWLINE_TABLE(comid);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${FLOWLINE_TABLE}_fcode_idx ON $FLOWLINE_TABLE(fcode);"

## Areas
areas=$(find . -iname nhdarea.shp)

set -- $areas
echo "Creating $AREA_TABLE schema"
shp2pgsql -p -D -t 2D -s $SRID $1 $AREA_TABLE | psql -d $DBNAME -q

for f in $areas; do
  echo "Importing $f..."
  shp2pgsql -a -D -t 2D -s $SRID -W $ENCODING $f $AREA_TABLE | pv | psql -d $DBNAME -q
done

psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${AREA_TABLE}_geom_gist ON $AREA_TABLE USING GIST(geom);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${AREA_TABLE}_ftype_idx ON $AREA_TABLE(ftype);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${AREA_TABLE}_comid_idx ON $AREA_TABLE(comid);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${AREA_TABLE}_fcode_idx ON $AREA_TABLE(fcode);"

## Waterbodies
waterbodies=$(find . -iname nhdwaterbody.shp)

set -- $waterbodies
echo "Creating $WB_TABLE schema"
shp2pgsql -p -D -t 2D -s $SRID $1 $WB_TABLE | psql -d $DBNAME -q

for f in $waterbodies; do
  echo "Importing $f..."
  shp2pgsql -a -D -t 2D -s $SRID -W $ENCODING $f $WB_TABLE | pv | psql -d $DBNAME -q
done

psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${WB_TABLE}_geom_gist ON $WB_TABLE USING GIST(geom);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${WB_TABLE}_ftype_idx ON $WB_TABLE(ftype);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${WB_TABLE}_comid_idx ON $WB_TABLE(comid);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${WB_TABLE}_fcode_idx ON $WB_TABLE(fcode);"

## Value-Added Attributes
vaas=$(find . -iname plusflowlinevaa.dbf)

set -- $vaas
echo "Creating plusflowlinevaa schema"
pgdbf -s LATIN1 $1 | psql -d $DBNAME -q
psql -d $DBNAME -c "TRUNCATE TABLE plusflowlinevaa;"

for f in $vaas; do
  echo "Importing $f..."
  pgdbf -CDs LATIN1 $f | pv | psql -d $DBNAME -q
done

psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY plusflowlinevaa_comid_idx ON plusflowlinevaa(comid);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY plusflowlinevaa_streamorde_idx ON plusflowlinevaa(streamorde);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY plusflowlinevaa_levelpathi_idx ON plusflowlinevaa(levelpathi);"

## fcode lookups

pgdbf -s LATIN1 $(find . -name NHDFcode.dbf | head -1) | pv | psql -d nhdplus -q
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY nhdfcode_fcode_idx ON nhdfcode(fcode);"

## Points
points=$(find . -iname nhdpoint.shp)

set -- $points
echo "Creating $PT_TABLE schema"
shp2pgsql -p -D -s $SRID $1 $PT_TABLE | psql -d $DBNAME -q

for f in $points; do
  echo "Importing $f..."
  shp2pgsql -a -D -s $SRID -W $ENCODING $f $PT_TABLE | pv | psql -d $DBNAME -q
done

psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${PT_TABLE}_geom_gist ON $PT_TABLE USING GIST(geom);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${PT_TABLE}_ftype_idx ON $PT_TABLE(ftype);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${PT_TABLE}_comid_idx ON $PT_TABLE(comid);"

## Lines
lines=$(find . -iname nhdline.shp)

set -- $lines
echo "Creating $LINE_TABLE schema"
shp2pgsql -p -D -s $SRID $1 $LINE_TABLE | psql -d $DBNAME -q

for f in $lines; do
  echo "Importing $f..."
  shp2pgsql -a -D -s $SRID -W $ENCODING $f $LINE_TABLE | pv | psql -d $DBNAME -q
done

psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${LINE_TABLE}_geom_gist ON $LINE_TABLE USING GIST(geom);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${LINE_TABLE}_ftype_idx ON $LINE_TABLE(ftype);"
psql -d $DBNAME -c "CREATE INDEX CONCURRENTLY ${LINE_TABLE}_comid_idx ON $LINE_TABLE(comid);"

# create table flowline_labels AS (SELECT levelpathi AS levelpathid, terminalpa AS terminalpathid, gnis_name AS name, streamorde AS streamorder, ST_SetSRID(ST_Simplify(ST_LineMerge(ST_Union(geom)), 0.01), 4269) AS geom FROM flowline LEFT JOIN plusflowlinevaa vaa ON vaa.comid=flowline.comid WHERE flowline.fcode != 46003 AND (streamorde IS NULL OR streamorde >= 3) AND gnis_name != '' GROUP BY levelpathi, terminalpa, gnis_name, streamorde, ftype);
# CREATE INDEX CONCURRENTLY flowline_labels_geom_gist ON flowline_labels USING GIST(geom);