source: subversion/applications/rendering/toposm/NED.py @ 27489

Last change on this file since 27489 was 25435, checked in by ahlzen, 9 years ago

Minor refactoring and cleanup of the TopOSM code.

File size: 5.6 KB
Line 
1#!/usr/bin/python
2
3"""NED.py: NED (National Elevation Dataset) preprocessing for TopOSM."""
4
5import numpy
6import mapnik
7from os import path, system
8from math import floor, ceil
9
10from env import *
11from common import *
12
13__author__      = "Lars Ahlzen"
14__copyright__   = "(c) Lars Ahlzen 2008-2011"
15__license__     = "GPLv2"
16
17# Size (side length, in degrees) of NED tiles that are cut up into
18# more manageable chunks.
19STEP = 0.5
20
21def getTilepath(basename):
22    return os.path.join(
23        NED13DIR, basename, 'grd' + basename + '_13', 'w001001.adf')
24
25def getTiles(envLL):
26    """Gets the (basename, envelope) of all (existing) 1/3 NED tiles
27    to cover the specified envelope"""
28    fromx = int(floor(envLL.minx))
29    tox = int(floor(envLL.maxx))
30    fromy = int(ceil(envLL.miny))
31    toy = int(ceil(envLL.maxy))
32    tiles = []
33    for y in range(fromy, toy+1):
34        for x in range(fromx, tox+1):
35            basename = 'n%02dw%03d' % (y, -x)
36            tilepath = getTilepath(basename)
37            if path.isfile(tilepath):
38                tiles.append((basename, mapnik.Envelope(x, y-1, x+1, y)))
39    return tiles
40
41def getSlice(prefix, x, y, suffix = '.tif'):
42    filename = prefix + '_%.1f_%.1f%s' % (float(x), float(y), suffix)
43    return path.join(TEMPDIR, filename)
44
45def getSlices(prefix, envLL, suffix = '.tif'):
46    fromx = floor(envLL.minx/STEP)*STEP
47    fromy = floor(envLL.miny/STEP)*STEP
48    tox = ceil(envLL.maxx/STEP)*STEP
49    toy = ceil(envLL.maxy/STEP)*STEP
50    slices = []
51    for y in numpy.arange(fromy, toy, STEP):
52        for x in numpy.arange(fromx, tox, STEP):
53            slicefile = getSlice(prefix, x, y, suffix)
54            if path.isfile(slicefile):
55                slices.append(slicefile)
56    return slices
57
58def getTilecoords(lat, lon, step = 1):
59    return (int(ceil(lat/float(step)))*float(step), \
60        int(floor(lon/float(step)))*float(step))
61
62def getTilename(lat, lon, step = 1):
63    (y, x) = get_ned13_tilecoords(lat, lon, step)
64    if step == 1:
65        return 'n%02dw%03d' % (y, -x)
66    return 'n%02.2fw%03.2f' % (y, -x)
67                       
68def prepDataFile(basename, env):
69    print 'Preparing NED 1/3" data file:', basename
70    print '  Converting to GeoTIFF...'
71    sourcepath = getTilepath(basename)
72    geotiff = path.join(TEMPDIR, basename + '.tif')
73    if not path.isfile(geotiff):
74        cmd = 'gdal_translate "%s" "%s"' % (sourcepath, geotiff)
75        #call(cmd, shell = True)
76        os.system(cmd)
77
78    print '  Generating contour lines...'
79    # split the GeoTIFF, since it's often too large otherwise
80    for y in numpy.arange(env.miny, env.maxy, STEP):
81        for x in numpy.arange(env.minx, env.maxx, STEP):
82            print '  Cutting geotiff slice...'
83            nedslice = getSlice('ned', x, y)
84            if not path.isfile(nedslice):
85                cmd = 'gdalwarp -q -te %f %f %f %f "%s" "%s"' % \
86                    (x, y, x+STEP, y+STEP, geotiff, nedslice)
87                #call(cmd, shell=True)
88                os.system(cmd)
89           
90            print '  Generating contour lines...'
91            contourbasefile = path.join(TEMPDIR, 'contours_' + str(x) + '_' + str(y))
92            contourfile = contourbasefile + '.shp'
93            if not path.isfile(contourfile):
94                cmd = 'gdal_contour -i %f -snodata 32767 -a height "%s" "%s"' % \
95                    (CONTOUR_INTERVAL, nedslice, contourfile)
96                #call(cmd, shell=True)
97                os.system(cmd)
98
99                print '  Importing contour lines...'
100                # NOTE: this assumes that the table is already set up
101                cmd = 'shp2pgsql -a -g way "%s" "%s" | psql -q "%s"' % \
102                    (contourfile, CONTOURS_TABLE, DATABASE)
103                #call(cmd, shell=True)
104                os.system(cmd)
105               
106                # Clear contents (but keep file to prevent us from importing
107                # these contours again).
108                writeEmpty(contourfile)
109               
110                # Remove shape index and attribute files.
111                tryRemove(contourbasefile + ".shx")
112                tryRemove(contourbasefile + ".dbf")
113
114            print '  Generating hillshade slice...'
115            hillshadeslice = getSlice('hillshade', x, y)
116            hillshadeslicePng = getSlice('hillshade', x, y, '.png')
117            if not path.isfile(hillshadeslicePng):
118                cmd = '"%s" "%s" "%s" -z 0.00001' % \
119                    (HILLSHADE, nedslice, hillshadeslice)
120                #call(cmd, shell = True)
121                os.system(cmd)
122                # convert to PNG + world file to save space
123                cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
124                    (hillshadeslice, hillshadeslicePng)
125                #call(cmd, shell = True)
126                os.system(cmd)
127                tryRemove(hillshadeslice)
128
129            print '  Generating colormap slice...'
130            colormapslice = getSlice('colormap', x, y)
131            colormapslicePng = getSlice('colormap', x, y, '.png')
132            if not path.isfile(colormapslicePng):
133                cmd = '"%s" "%s" "%s" "%s"' % \
134                    (COLORRELIEF, nedslice, COLORFILE, colormapslice)
135                #call(cmd, shell = True)
136                os.system(cmd)
137                # convert to PNG + world file to save space
138                cmd = 'gdal_translate -quiet -of PNG -co WORLDFILE=YES "%s" "%s"' % \
139                    (colormapslice, colormapslicePng)
140                #call(cmd, shell = True)
141                os.system(cmd)
142                tryRemove(colormapslice)
143
144            writeEmpty(nedslice) # don't need the raw slice anymore.
145               
146    writeEmpty(geotiff) # done with this GeoTIFF.
147
Note: See TracBrowser for help on using the repository browser.