Added NLCD preprocessing.

Improved contour data preprocessing.

NED.py

@@ -72,19 +72,34 @@ def getTilename(lat, lon, step = 1):
         return 'n%02dw%03d' % (y, -x)
     return 'n%02.2fw%03.2f' % (y, -x) 
 
-def simplifyContours():
-    cmd = 'UPDATE %s SET way = ST_Simplify(way, 1.0);' % (CONTOURS_TABLE)
-    # TODO: Finish
+def convertContourElevationsToFt():
+
+    templ = 'UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT)' + \
+        ' WHERE height_ft IS NULL'
+    sql = templ % (CONTOURS_TABLE)
+    runSql(sql)
+
+# NOTE: This removes a LOT of artifacts around shorelines. Unfortunately,
+# it also removes any legitimate 0ft contours (which may exist around
+# areas below sea-level).
+# TODO: Find a better workaround.
+def removeSeaLevelContours():
+    sql = 'DELETE FROM %s WHERE height_ft = 0;' % (CONTOURS_TABLE)
+    runSql(sql)
+
+def simplifyContours(level = 1.0):
+    sql = 'UPDATE %s SET way = ST_Simplify(way, %f);' % (CONTOURS_TABLE, level)
+    runSql(sql)
+
+def clusterContoursOnGeoColumn():
+    sql = 'CLUSTER %s_way_gist ON %s;' % (CONTOURS_TABLE, CONTOURS_TABLE)
+    runSql(sql)
+
+def analyzeContoursTable():
+    sql = 'ANALYZE %s;' % (CONTOURS_TABLE)
+    runSql(sql)
 
 def prepDataFile(basename, env):
-    #print 'Preparing NED 1/3" data file:', basename
-    #print '  Converting to GeoTIFF...'
-    #sourcepath = getTilepath(basename)
-    #geotiff = path.join(TEMPDIR, basename + '.tif')
-    #if not path.isfile(geotiff):
-    #    cmd = 'gdal_translate "%s" "%s"' % (sourcepath, geotiff)
-    #    #call(cmd, shell = True)
-    #    os.system(cmd)
     geotiff = getTilepath(basename)
 
     # split the GeoTIFF, since it's often too large otherwise
@@ -126,33 +141,5 @@ def prepDataFile(basename, env):
                 tryRemove(contourbasefile + ".shx")
                 tryRemove(contourbasefile + ".dbf")
 
-            #hillshadeslice = getSlice('hillshade', x, y)
-            #hillshadeslicePng = getSlice('hillshade', x, y, '.png')
-            #if not path.isfile(hillshadeslicePng):
-            #    print '  Generating hillshade slice...'
-            #    cmd = '"%s" "%s" "%s" -z 0.00001' % \
-            #        (HILLSHADE, nedslice, hillshadeslice)
-            #    os.system(cmd)
-            #    # convert to PNG + world file to save space
-            #    cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
-            #        (hillshadeslice, hillshadeslicePng)
-            #    os.system(cmd)
-            #    tryRemove(hillshadeslice)
-
-            #colormapslice = getSlice('colormap', x, y)
-            #colormapslicePng = getSlice('colormap', x, y, '.png')
-            #if not path.isfile(colormapslicePng):
-            #    print '  Generating colormap slice...'
-            #    cmd = '"%s" "%s" "%s" "%s"' % \
-            #        (COLORRELIEF, nedslice, COLORFILE, colormapslice)
-            #    os.system(cmd)
-            #    # convert to PNG + world file to save space
-            #    cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
-            #        (colormapslice, colormapslicePng)
-            #    os.system(cmd)
-            #    tryRemove(colormapslice)
-
             writeEmpty(nedslice) # don't need the raw slice anymore.
-
-    #writeEmpty(geotiff) # done with this GeoTIFF.
-
+

common.py

@@ -69,4 +69,8 @@ def writeEmpty(filename):
     finally:
         fslock.release()
 
+def runSql(sql):
+    command = "psql -d %s -q -c \"%s\"" % (DATABASE, sql)
+    print command
+    os.system(command)
 

prep_contours_table

@@ -9,4 +9,3 @@ DBCMD="psql -q -p $DB_PORT $DB_NAME $DB_USER"
 
 shp2pgsql -p -I -g way contours.shp $CONTOURS_TABLE | $DBCMD
 echo "ALTER TABLE $CONTOURS_TABLE ADD COLUMN height_ft INT" | $DBCMD
-

prep_toposm_data

@@ -83,3 +83,13 @@ for d in $COLORMAP_DIR $HILLSHADE_DIR $HYPSORELIEF_DIR ; do
     gdalbuildvrt all.vrt *.tif
     popd
 done
+
+# reproject NLCD data, add overviews, and copy VRT file (with
+# the custom color table).
+gdalwarp -multi-of GTiff \
+    -co ZLEVEL=9 -co COMPRESS=DEFLATE -co PREDICTOR=2 -co BIGTIFF=YES \
+    -t_srs EPSG:900913 $NLCD_DIR/nlcd2006_landcover_2-14-11.img \
+    $NLCD_DIR/nlcd2006.tif
+gdaladdo $NLCD_DIR/nlcd2006.tif 2 4 8 16 32 64
+cp nlcd2006.vrt $NLCD_DIR
+

set-toposm-env

@@ -14,6 +14,7 @@ export NHD_TABLE_PREFIX="nhd"
 # Data directories
 export WORLD_BOUNDARIES_DIR="geodata/osm/world_boundaries"
 export NHD_DIR="/home/lars/titan/projects/toposm/geodata/NHD" # NHD shapefiles
+export NLCD_DIR="geodata/nlcd2006"
 export NED13_DIR="geodata/ned13" # NED 1/3" GeoTIFFs
 export HILLSHADE_DIR="geodata/hillshade" # Hillshade GeoTIFFs
 export COLORMAP_DIR="geodata/colormap" # Colormap GeoTIFFs

toposm.py

@@ -244,20 +244,12 @@ def prepareData(envLLs):
             manager.addJob("Preparing %s" % (tile[0]), NED.prepDataFile, tile)
     manager.finish()
 
-    console.printMessage("Converting m to ft")
-    templ = 'echo "UPDATE %s SET height_ft = CAST(height * 3.28085 AS INT) ' \
-        + 'WHERE height_ft IS NULL;" | psql -q "%s"'
-    cmd = templ % (CONTOURS_TABLE, DATABASE)
-    os.system(cmd)
-
-    # NOTE: This removes a LOT of artifacts around shorelines. Unfortunately,
-    # it also removes any legitimate 0ft contours (which may exist around
-    # areas below sea-level).
-    # TODO: Find a better workaround.
-    console.printMessage('Removing sea-level contours')
-    templ = 'echo "DELETE FROM %s WHERE height_ft = 0;" | psql -q "%s"'
-    cmd = templ % (CONTOURS_TABLE, DATABASE)
-    os.system(cmd)
+    console.printMessage("Postprocessing contours...")
+    NED.removeSeaLevelContours()
+    NED.simplifyContours(1.0)
+    NED.convertContourElevationsToFt()
+    NED.clusterContoursOnGeoColumn()
+    NED.analyzeContoursTable()
 
 def renderTiles(envLLs, minz, maxz):
     if not hasattr(envLLs, '__iter__'):
@@ -320,26 +312,22 @@ def printSyntax():
     print "Areas are named entities in areas.py."
 
 if __name__ == "__main__":
-    try:
-        cmd = sys.argv[1]
-        if cmd == 'render':
-            areaname = sys.argv[2]
-            minzoom = int(sys.argv[3])
-            maxzoom = int(sys.argv[4])
-            env = vars(areas)[areaname]            
-            print "Render: %s %s, z: %d-%d" % (areaname, env, minzoom, maxzoom)
-            BASE_TILE_DIR = path.join(BASE_TILE_DIR, areaname)
-            renderTiles(env, minzoom, maxzoom)
-        elif cmd == 'prep':
-            areaname = sys.argv[2]
-            env = vars(areas)[areaname]
-            print "Prepare data: %s %s" % (areaname, env)
-            prepareData(env)
-        elif cmd == 'info':
-            toposmInfo()
-        else:
-            raise Exception()
-    except:
+    cmd = sys.argv[1]
+    if cmd == 'render':
+        areaname = sys.argv[2]
+        minzoom = int(sys.argv[3])
+        maxzoom = int(sys.argv[4])
+        env = vars(areas)[areaname]            
+        print "Render: %s %s, z: %d-%d" % (areaname, env, minzoom, maxzoom)
+        BASE_TILE_DIR = path.join(BASE_TILE_DIR, areaname)
+        renderTiles(env, minzoom, maxzoom)
+    elif cmd == 'prep':
+        areaname = sys.argv[2]
+        env = vars(areas)[areaname]
+        print "Prepare data: %s %s" % (areaname, env)
+        prepareData(env)
+    elif cmd == 'info':
+        toposmInfo()
+    else:
         printSyntax()
         sys.exit(1)
-