| 1 | #!/usr/bin/python
|
|---|
| 2 |
|
|---|
| 3 | # Lakewalker II.
|
|---|
| 4 | # 2007-07-09: Initial public release (v0.2) by Dshpak
|
|---|
| 5 | # 2007-07-10: v0.3: Added support for OSM input files. Also reduced start_radius_big, and fixed a division-by-zero error in the degenerate case of point_line_distance().
|
|---|
| 6 | # 2007-08-04: Added experimental non-recursive Douglas-Peucker algorithm (--dp-nr). Added --startdir option.
|
|---|
| 7 | # 2007-08-06: v0.4: Bounding box support (--left, --right, --top, and --bottom), as well as new --josm mode for JOSM integration. This isn't perfect yet, but it's a good start.
|
|---|
| 8 |
|
|---|
| 9 | """Lakewalker II - A tool to automatically download and trace Landsat imagery, to generate OpenStreetMap data.
|
|---|
| 10 |
|
|---|
| 11 | Requires the Python Imaging Library, available from http://www.pythonware.com/products/pil/"""
|
|---|
| 12 |
|
|---|
| 13 | version="Lakewalker II v0.4"
|
|---|
| 14 |
|
|---|
| 15 | # TODO:
|
|---|
| 16 | # - Accept threshold tags in OSM input files.
|
|---|
| 17 | # - Command line options:
|
|---|
| 18 | # - direction (deosil/widdershins)
|
|---|
| 19 | # - layer to use (monochrome only?)
|
|---|
| 20 | # - additional tags for the way
|
|---|
| 21 | # - Landsat download retries (count and delay)
|
|---|
| 22 | # - radii for loop detection
|
|---|
| 23 | # - fname for z12 tile output (list or script)
|
|---|
| 24 | # - path to tilesGen for z12 script output
|
|---|
| 25 | # - debug mode (extra tags on nodes) (make this non-uploadable?)
|
|---|
| 26 | # - The "got stuck" message should output coords
|
|---|
| 27 | # - Better "got stuck" detection (detect mini loops)
|
|---|
| 28 | # - Offset nodes outwards by a half-pixel or so, to prevent duplicate segments
|
|---|
| 29 | # - Automatic threshold detection
|
|---|
| 30 | # - (Correct/tested) non-recursive DP simplification
|
|---|
| 31 | #
|
|---|
| 32 | #
|
|---|
| 33 | # For JOSM integration:
|
|---|
| 34 | # - Add a --tmpdir option that controls the location of the Landsat
|
|---|
| 35 | # tiles, the results text file, and the lake.osm file (unless it's
|
|---|
| 36 | # got an absolute path
|
|---|
| 37 | # - Add a --nocache option that keeps WMS tiles in memory but not on disk.
|
|---|
| 38 |
|
|---|
| 39 | import math
|
|---|
| 40 | import os
|
|---|
| 41 | import urllib
|
|---|
| 42 | from PIL import Image
|
|---|
| 43 | import OSMReader
|
|---|
| 44 | import time
|
|---|
| 45 | import optparse
|
|---|
| 46 |
|
|---|
| 47 | options = None
|
|---|
| 48 |
|
|---|
| 49 | start_radius_big = 0.001
|
|---|
| 50 | start_radius_small = 0.0002
|
|---|
| 51 |
|
|---|
| 52 | dirs = ((1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1))
|
|---|
| 53 | dirnames = ["east", "northeast", "north", "northwest", "west", "southwest", "south", "southeast"]
|
|---|
| 54 | dirabbrevs = ["e", "ne", "n", "nw", "w", "sw", "s", "se"]
|
|---|
| 55 |
|
|---|
| 56 | class FatalError(Exception):
|
|---|
| 57 | pass
|
|---|
| 58 |
|
|---|
| 59 | def message(s):
|
|---|
| 60 | if options.josm_mode:
|
|---|
| 61 | print "m %s" % s
|
|---|
| 62 | else:
|
|---|
| 63 | print s
|
|---|
| 64 |
|
|---|
| 65 | def error(s):
|
|---|
| 66 | if options.josm_mode:
|
|---|
| 67 | print "e %s" % s
|
|---|
| 68 | else:
|
|---|
| 69 | print s
|
|---|
| 70 |
|
|---|
| 71 | class BBox:
|
|---|
| 72 | def __init__(self, top = 90, left = -180, bottom = -90, right = 180):
|
|---|
| 73 | self.left = left
|
|---|
| 74 | self.right = right
|
|---|
| 75 | self.top = top
|
|---|
| 76 | self.bottom = bottom
|
|---|
| 77 | def contains(self, loc):
|
|---|
| 78 | (lat, lon) = loc
|
|---|
| 79 | if lat > self.top or lat < self.bottom:
|
|---|
| 80 | return False
|
|---|
| 81 | if (self.right - self.left) % 360 == 0:
|
|---|
| 82 | return True
|
|---|
| 83 | return (lon - self.left) % 360 <= (self.right - self.left) % 360
|
|---|
| 84 |
|
|---|
| 85 | def download_landsat(c1, c2, width, height, fname):
|
|---|
| 86 | layer = "global_mosaic_base"
|
|---|
| 87 | style = options.wmslayer
|
|---|
| 88 |
|
|---|
| 89 | (min_lat, min_lon) = c1
|
|---|
| 90 | (max_lat, max_lon) = c2
|
|---|
| 91 |
|
|---|
| 92 | message("Downloading Landsat tile for (%.6f,%.6f)-(%.6f,%.6f)." % (min_lat, min_lon, max_lat, max_lon))
|
|---|
| 93 | url = "http://onearth.jpl.nasa.gov/wms.cgi?request=GetMap&layers=%s&styles=%s&srs=EPSG:4326&format=image/png&bbox=%0.6f,%0.6f,%0.6f,%0.6f&width=%d&height=%d" % (layer, style, min_lon, min_lat, max_lon, max_lat, width, height)
|
|---|
| 94 | #print url
|
|---|
| 95 | try:
|
|---|
| 96 | urllib.urlretrieve(url, fname)
|
|---|
| 97 | except IOError, e:
|
|---|
| 98 | raise FatalError("Error downloading tile: %s" % e.strerror)
|
|---|
| 99 | if not os.path.exists(fname):
|
|---|
| 100 | raise FatalError("Error: Could not retreive url %s" % url)
|
|---|
| 101 |
|
|---|
| 102 | def xy_to_geo(xy):
|
|---|
| 103 | (x, y) = xy
|
|---|
| 104 | (lat, lon) = (y / float(options.resolution), x / float(options.resolution))
|
|---|
| 105 | return (lat, lon)
|
|---|
| 106 |
|
|---|
| 107 | def geo_to_xy(geo):
|
|---|
| 108 | (lat, lon) = geo
|
|---|
| 109 | coord = lambda L: math.floor(L * options.resolution + 0.5)
|
|---|
| 110 | (x, y) = (coord(lon), coord(lat))
|
|---|
| 111 | return (x, y)
|
|---|
| 112 |
|
|---|
| 113 | class WMSManager:
|
|---|
| 114 | def __init__(self):
|
|---|
| 115 | self.images = {}
|
|---|
| 116 |
|
|---|
| 117 | def get_tile(self, xy):
|
|---|
| 118 | fail_count = 0
|
|---|
| 119 | im = None
|
|---|
| 120 | while im is None and fail_count < 4:
|
|---|
| 121 | (x, y) = xy
|
|---|
| 122 | bottom_left_xy = (int(math.floor(x / options.tilesize)) * options.tilesize, int(math.floor(y / options.tilesize)) * options.tilesize)
|
|---|
| 123 | top_right_xy = (bottom_left_xy[0] + options.tilesize, bottom_left_xy[1] + options.tilesize)
|
|---|
| 124 | fname = options.wmslayer+"/landsat_%d_%d_xy_%d_%d.png" % (options.resolution, options.tilesize, bottom_left_xy[0], bottom_left_xy[1])
|
|---|
| 125 | im = self.images.get(fname, None)
|
|---|
| 126 | if im is None:
|
|---|
| 127 | if not os.path.exists(fname):
|
|---|
| 128 | bottom_left = xy_to_geo(bottom_left_xy)
|
|---|
| 129 | top_right = xy_to_geo(top_right_xy)
|
|---|
| 130 | download_landsat(bottom_left, top_right, options.tilesize, options.tilesize, fname)
|
|---|
| 131 | if not os.path.exists(fname):
|
|---|
| 132 | raise FatalError("Error: Could not get image file %s" % fname)
|
|---|
| 133 | try:
|
|---|
| 134 | im = Image.open(fname)
|
|---|
| 135 | self.images[fname] = im
|
|---|
| 136 | message("Using imagery in %s for %s." % (fname, xy_to_geo(xy)))
|
|---|
| 137 | except IOError:
|
|---|
| 138 | error("Download was corrupt...Deleting %s..." % fname)
|
|---|
| 139 | os.unlink(fname)
|
|---|
| 140 | im = None
|
|---|
| 141 |
|
|---|
| 142 | if im is None:
|
|---|
| 143 | message("Sleeping and retrying download...")
|
|---|
| 144 | time.sleep(4)
|
|---|
| 145 | fail_count = fail_count + 1
|
|---|
| 146 |
|
|---|
| 147 | if im is None:
|
|---|
| 148 | #if os.path.exists(fname):
|
|---|
| 149 | # print open(fname).readlines()
|
|---|
| 150 | raise FatalError("Couldn't get image file %s." % fname)
|
|---|
| 151 |
|
|---|
| 152 | #return (im, top_left_xy)
|
|---|
| 153 | return (im, bottom_left_xy)
|
|---|
| 154 |
|
|---|
| 155 | def get_pixel(self, xy):
|
|---|
| 156 | (x, y) = xy
|
|---|
| 157 | (tile, (tx, ty)) = self.get_tile(xy)
|
|---|
| 158 | tile_xy = (x - tx, (options.tilesize - 1) - (y - ty))
|
|---|
| 159 | #print "%s maps to %s" % (xy, tile_xy)
|
|---|
| 160 | return tile.getpixel(tile_xy)
|
|---|
| 161 |
|
|---|
| 162 | def trace_lake(loc, threshold, start_dir, bbox):
|
|---|
| 163 | wms = WMSManager()
|
|---|
| 164 | xy = geo_to_xy(loc)
|
|---|
| 165 | nodelist = []
|
|---|
| 166 |
|
|---|
| 167 | message("Starting coordinate: %.4f, %.4f" % loc)
|
|---|
| 168 | message("Starting position: %d, %d" % xy)
|
|---|
| 169 |
|
|---|
| 170 | if not bbox.contains(loc):
|
|---|
| 171 | raise FatalError("Error: Starting location is outside bounding box!")
|
|---|
| 172 |
|
|---|
| 173 | while True:
|
|---|
| 174 | loc = xy_to_geo(xy)
|
|---|
| 175 | if not bbox.contains(loc):
|
|---|
| 176 | break
|
|---|
| 177 |
|
|---|
| 178 | v = wms.get_pixel(xy)
|
|---|
| 179 | if v > threshold:
|
|---|
| 180 | break
|
|---|
| 181 |
|
|---|
| 182 | xy = (xy[0] + dirs[start_dir][0], xy[1] + dirs[start_dir][1])
|
|---|
| 183 |
|
|---|
| 184 | start_xy = xy
|
|---|
| 185 | start_loc = xy_to_geo(xy)
|
|---|
| 186 | message("Found shore at lat %.4f lon %.4f" % start_loc)
|
|---|
| 187 |
|
|---|
| 188 | #dirs = ((1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1), (0, 1), (1, 1))
|
|---|
| 189 | last_dir = start_dir
|
|---|
| 190 |
|
|---|
| 191 | detect_loop = False
|
|---|
| 192 |
|
|---|
| 193 | for i in xrange(options.maxnodes):
|
|---|
| 194 | if i % 250 == 0:
|
|---|
| 195 | if i > 0:
|
|---|
| 196 | message("%s nodes so far..." % i)
|
|---|
| 197 |
|
|---|
| 198 | for d in xrange(1, len(dirs)):
|
|---|
| 199 | new_dir = dirs[(last_dir + d + 4) % 8]
|
|---|
| 200 | test_xy = (xy[0] + new_dir[0], xy[1] + new_dir[1])
|
|---|
| 201 | test_loc = xy_to_geo(test_xy)
|
|---|
| 202 | if not bbox.contains(test_loc):
|
|---|
| 203 | break
|
|---|
| 204 |
|
|---|
| 205 | v = wms.get_pixel(test_xy)
|
|---|
| 206 | #print "%s: %s: %s" % (new_dir, test_xy, v)
|
|---|
| 207 | if v > threshold:
|
|---|
| 208 | break
|
|---|
| 209 |
|
|---|
| 210 | if d == 8:
|
|---|
| 211 | error("Got stuck.")
|
|---|
| 212 | break
|
|---|
| 213 |
|
|---|
| 214 | #print "Moving to %s, direction %s (was %s)" % (test_xy, new_dir, dirs[last_dir])
|
|---|
| 215 | last_dir = (last_dir + d + 4) % 8
|
|---|
| 216 | xy = test_xy
|
|---|
| 217 |
|
|---|
| 218 | if xy == start_xy:
|
|---|
| 219 | break
|
|---|
| 220 |
|
|---|
| 221 | loc = xy_to_geo(xy)
|
|---|
| 222 | nodelist.append(loc)
|
|---|
| 223 |
|
|---|
| 224 | start_proximity = (loc[0] - start_loc[0]) ** 2 + (loc[1] - start_loc[1]) ** 2
|
|---|
| 225 | if detect_loop:
|
|---|
| 226 | if start_proximity < start_radius_small ** 2:
|
|---|
| 227 | break
|
|---|
| 228 | else:
|
|---|
| 229 | if start_proximity > start_radius_big ** 2:
|
|---|
| 230 | detect_loop = True
|
|---|
| 231 | return nodelist
|
|---|
| 232 |
|
|---|
| 233 | def vertex_reduce(nodes, proximity):
|
|---|
| 234 | test_v = nodes[0]
|
|---|
| 235 | reduced_nodes = [test_v]
|
|---|
| 236 | prox_sq = proximity ** 2
|
|---|
| 237 | for v in nodes:
|
|---|
| 238 | if (v[0] - test_v[0]) ** 2 + (v[1] - test_v[1]) ** 2 > prox_sq:
|
|---|
| 239 | reduced_nodes.append(v)
|
|---|
| 240 | test_v = v
|
|---|
| 241 | return reduced_nodes
|
|---|
| 242 |
|
|---|
| 243 | def point_line_distance(p0, p1, p2):
|
|---|
| 244 | ((x0, y0), (x1, y1), (x2, y2)) = (p0, p1, p2)
|
|---|
| 245 |
|
|---|
| 246 | if x2 == x1 and y2 == y1:
|
|---|
| 247 | # Degenerate cast: the "line" is actually a point.
|
|---|
| 248 | return math.sqrt((x1-x0)**2 + (y1-y0)**2)
|
|---|
| 249 | else:
|
|---|
| 250 | # I don't understand this at all. Thank you, Mathworld.
|
|---|
| 251 | # http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
|
|---|
| 252 | return abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / math.sqrt((x2-x1)**2 + (y2-y1)**2)
|
|---|
| 253 |
|
|---|
| 254 | def douglas_peucker(nodes, epsilon):
|
|---|
| 255 | #print "Running DP on %d nodes" % len(nodes)
|
|---|
| 256 | farthest_node = None
|
|---|
| 257 | farthest_dist = 0
|
|---|
| 258 | first = nodes[0]
|
|---|
| 259 | last = nodes[-1]
|
|---|
| 260 |
|
|---|
| 261 | for i in xrange(1, len(nodes) - 1):
|
|---|
| 262 | d = point_line_distance(nodes[i], first, last)
|
|---|
| 263 | if d > farthest_dist:
|
|---|
| 264 | farthest_dist = d
|
|---|
| 265 | farthest_node = i
|
|---|
| 266 |
|
|---|
| 267 | if farthest_dist > epsilon:
|
|---|
| 268 | seg_a = douglas_peucker(nodes[0:farthest_node+1], epsilon)
|
|---|
| 269 | seg_b = douglas_peucker(nodes[farthest_node:-1], epsilon)
|
|---|
| 270 | #print "Minimized %d nodes to %d + %d nodes" % (len(nodes), len(seg_a), len(seg_b))
|
|---|
| 271 | nodes = seg_a[:-1] + seg_b
|
|---|
| 272 | else:
|
|---|
| 273 | return [nodes[0], nodes[-1]]
|
|---|
| 274 |
|
|---|
| 275 | return nodes
|
|---|
| 276 |
|
|---|
| 277 | def dp_findpoint(nodes, start, end):
|
|---|
| 278 | farthest_node = None
|
|---|
| 279 | farthest_dist = 0
|
|---|
| 280 | #print "dp_findpoint(nodes, %s, %s)" % (start, end)
|
|---|
| 281 | first = nodes[start]
|
|---|
| 282 | last = nodes[end]
|
|---|
| 283 |
|
|---|
| 284 | for i in xrange(start + 1, end):
|
|---|
| 285 | d = point_line_distance(nodes[i], first, last)
|
|---|
| 286 | if d > farthest_dist:
|
|---|
| 287 | farthest_dist = d
|
|---|
| 288 | farthest_node = i
|
|---|
| 289 |
|
|---|
| 290 | return (farthest_node, farthest_dist)
|
|---|
| 291 |
|
|---|
| 292 | def douglas_peucker_nonrecursive(nodes, epsilon):
|
|---|
| 293 | #print "Running DP on %d nodes" % len(nodes)
|
|---|
| 294 | command_stack = [(0, len(nodes) - 1)]
|
|---|
| 295 | result_stack = []
|
|---|
| 296 |
|
|---|
| 297 | while len(command_stack) > 0:
|
|---|
| 298 | cmd = command_stack.pop()
|
|---|
| 299 | if type(cmd) == tuple:
|
|---|
| 300 | (start, end) = cmd
|
|---|
| 301 | (node, dist) = dp_findpoint(nodes, start, end)
|
|---|
| 302 | if dist > epsilon:
|
|---|
| 303 | command_stack.append("+")
|
|---|
| 304 | command_stack.append((start, node))
|
|---|
| 305 | command_stack.append((node, end))
|
|---|
| 306 | else:
|
|---|
| 307 | result_stack.append((start, end))
|
|---|
| 308 | elif cmd == "+":
|
|---|
| 309 | first = result_stack.pop()
|
|---|
| 310 | second = result_stack.pop()
|
|---|
| 311 | if first[-1] == second[0]:
|
|---|
| 312 | result_stack.append(first + second[1:])
|
|---|
| 313 | #print "Added %s and %s; result is %s" % (first, second, result_stack[-1])
|
|---|
| 314 | else:
|
|---|
| 315 | error("ERROR: Cannot connect nodestrings!")
|
|---|
| 316 | #print first
|
|---|
| 317 | #print second
|
|---|
| 318 | return
|
|---|
| 319 | else:
|
|---|
| 320 | error("ERROR: Can't understand command \"%s\"" % (cmd,))
|
|---|
| 321 | return
|
|---|
| 322 |
|
|---|
| 323 | if len(result_stack) == 1:
|
|---|
| 324 | return [nodes[x] for x in result_stack[0]]
|
|---|
| 325 | else:
|
|---|
| 326 | error("ERROR: Command stack is empty but result stack has %d nodes!" % len(result_stack))
|
|---|
| 327 | return
|
|---|
| 328 |
|
|---|
| 329 | farthest_node = None
|
|---|
| 330 | farthest_dist = 0
|
|---|
| 331 | first = nodes[0]
|
|---|
| 332 | last = nodes[-1]
|
|---|
| 333 |
|
|---|
| 334 | for i in xrange(1, len(nodes) - 1):
|
|---|
| 335 | d = point_line_distance(nodes[i], first, last)
|
|---|
| 336 | if d > farthest_dist:
|
|---|
| 337 | farthest_dist = d
|
|---|
| 338 | farthest_node = i
|
|---|
| 339 |
|
|---|
| 340 | if farthest_dist > epsilon:
|
|---|
| 341 | seg_a = douglas_peucker(nodes[0:farthest_node+1], epsilon)
|
|---|
| 342 | seg_b = douglas_peucker(nodes[farthest_node:-1], epsilon)
|
|---|
| 343 | #print "Minimized %d nodes to %d + %d nodes" % (len(nodes), len(seg_a), len(seg_b))
|
|---|
| 344 | nodes = seg_a[:-1] + seg_b
|
|---|
| 345 | else:
|
|---|
| 346 | return [nodes[0], nodes[-1]]
|
|---|
| 347 |
|
|---|
| 348 | return nodes
|
|---|
| 349 |
|
|---|
| 350 | def output_to_josm(lakelist):
|
|---|
| 351 | # Description of JOSM output format:
|
|---|
| 352 | # m text - Status message text, to be displayed in a display window
|
|---|
| 353 | # e text - Error message text
|
|---|
| 354 | # s nnn - Start full node list, nnn tracings following
|
|---|
| 355 | # t nnn - Start tracing node list, nnn nodes following (i.e. there will be one of these for each lake where multiple start points specified)
|
|---|
| 356 | # n lat lon - A node
|
|---|
| 357 | # x - End of data
|
|---|
| 358 | print "s %s" % len(lakelist)
|
|---|
| 359 | for nodelist in lakelist:
|
|---|
| 360 | print "t %s" % len(nodelist)
|
|---|
| 361 | for node in nodelist:
|
|---|
| 362 | print "n %.7f %.7f" % (node[0], node[1])
|
|---|
| 363 | print "x"
|
|---|
| 364 |
|
|---|
| 365 | def write_osm(f, lakelist, waysize):
|
|---|
| 366 | f.write('<osm version="0.4">')
|
|---|
| 367 | cur_id = -1
|
|---|
| 368 | way_count = 0
|
|---|
| 369 | for nodelist in lakelist:
|
|---|
| 370 | first_node_id = cur_id
|
|---|
| 371 | cur_way = []
|
|---|
| 372 | ways = [cur_way]
|
|---|
| 373 | for loc in nodelist:
|
|---|
| 374 | #f.write(' <node id="%d" lat="%.6f" lon="%.6f"><tag k="order" v="%d"/><tag k="x" v="%d"/><tag k="y" v="%d"/></node>\n' % (cur_id, loc[0], loc[1], -cur_id, geo_to_xy(loc)[0], geo_to_xy(loc)[1]))
|
|---|
| 375 | f.write('<node id="%d" lat="%.6f" lon="%.6f"/>' % (cur_id, loc[0], loc[1]))
|
|---|
| 376 | last_node_id = cur_id
|
|---|
| 377 | cur_id = cur_id - 1
|
|---|
| 378 |
|
|---|
| 379 | # print "Nodes: %d, %d" % (first_node_id, last_node_id)
|
|---|
| 380 |
|
|---|
| 381 | first_segment_id = cur_id
|
|---|
| 382 | for seg in xrange(first_node_id, last_node_id, -1):
|
|---|
| 383 | f.write('<segment id="%d" from="%d" to="%d"/>' % (cur_id, seg, seg-1))
|
|---|
| 384 | cur_id = cur_id - 1
|
|---|
| 385 | f.write('<segment id="%d" from="%d" to="%d"/>' % (cur_id, last_node_id, first_node_id))
|
|---|
| 386 | last_segment_id = cur_id
|
|---|
| 387 | cur_id = cur_id - 1
|
|---|
| 388 |
|
|---|
| 389 | # print "Segments: %d, %d" % (first_segment_id, last_segment_id)
|
|---|
| 390 |
|
|---|
| 391 | for seg in xrange(first_segment_id, last_segment_id - 1, -1):
|
|---|
| 392 | if len(cur_way) >= waysize:
|
|---|
| 393 | cur_way = []
|
|---|
| 394 | ways.append(cur_way)
|
|---|
| 395 | cur_way.append(seg)
|
|---|
| 396 | for way in ways:
|
|---|
| 397 | f.write('<way id="%d">' % cur_id)
|
|---|
| 398 | cur_id = cur_id - 1
|
|---|
| 399 | for seg in way:
|
|---|
| 400 | f.write('<seg id="%d"/>' % seg)
|
|---|
| 401 | f.write('<tag k="natural" v="%s"/>' % options.natural_type)
|
|---|
| 402 | f.write('<tag k="source" v="Dshpak_landsat_lakes"/>')
|
|---|
| 403 | f.write('</way>')
|
|---|
| 404 | way_count = way_count + len(ways)
|
|---|
| 405 |
|
|---|
| 406 | f.write('</osm>')
|
|---|
| 407 | message("Generated %d %s." % (way_count, ["way", "ways"][way_count > 0]))
|
|---|
| 408 |
|
|---|
| 409 | def get_locs(infile):
|
|---|
| 410 | nodes = []
|
|---|
| 411 | segments = []
|
|---|
| 412 | ways = []
|
|---|
| 413 | reader = OSMReader.OSMReader()
|
|---|
| 414 | reader.nodeHandler = lambda x: nodes.append(x)
|
|---|
| 415 | reader.segmentHandler = lambda x: segments.append(x)
|
|---|
| 416 | reader.wayHandler = lambda x: ways.append(x)
|
|---|
| 417 | reader.run(file(infile))
|
|---|
| 418 | if len(segments) > 0 or len(ways) > 0:
|
|---|
| 419 | raise FatalError("Error: Input file must only contain nodes -- no segments or ways.")
|
|---|
| 420 | if len(nodes) == 0:
|
|---|
| 421 | raise FatalError("Error: No nodes found in input file.")
|
|---|
| 422 | return [((node.lat, node.lon), options.threshold) for node in nodes]
|
|---|
| 423 |
|
|---|
| 424 | def main():
|
|---|
| 425 | parser = optparse.OptionParser(version=version)
|
|---|
| 426 |
|
|---|
| 427 | parser.add_option("--lat", type="float", metavar="LATITUDE", help="Starting latitude. Required, unless --infile is used..")
|
|---|
| 428 | parser.add_option("--lon", type="float", metavar="LONGITUDE", help="Starting longitude. Required, unless --infile is used.")
|
|---|
| 429 | parser.add_option("--startdir", type="string", default="east", metavar="DIR", help="Direction to travel from start position when seeking land. Defaults to \"east\".")
|
|---|
| 430 | parser.add_option("--infile", "-i", type="string", metavar="FILE", help="OSM file containing nodes representing starting points.")
|
|---|
| 431 | parser.add_option("--out", "-o", default="lake.osm", dest="outfile", metavar="FILE", help="Output filename. Defaults to lake.osm.")
|
|---|
| 432 | parser.add_option("--threshold", "-t", type="int", default="35", metavar="VALUE", help="Maximum gray value to accept as water (based on Landsat IR-1 data). Can be in the range 0-255. Defaults to 35.")
|
|---|
| 433 | parser.add_option("--maxnodes", type="int", default="50000", metavar="N", help="Maximum number of nodes to generate before bailing out. Defaults to 50000.")
|
|---|
| 434 | parser.add_option("--waylength", type="int", default=250, metavar="MAXLEN", help="Maximum nuber of nodes allowed in one way. Defaults to 250.")
|
|---|
| 435 | parser.add_option("--landsat-res", type="int", default=4000, dest="resolution", metavar="RES", help="Resolution of Landsat tiles, measured in pixels per degree. Defaults to 4000.")
|
|---|
| 436 | parser.add_option("--tilesize", type="int", default=2000, help="Size of one landsat tile, measured in pixels. Defaults to 2000.")
|
|---|
| 437 | parser.add_option("--no-dp", action="store_false", dest="use_dp", default=True, help="Disable Douglas-Peucker line simplification (not recommended)")
|
|---|
| 438 | parser.add_option("--dp-epsilon", type="float", metavar="EPSILON", default=0.0003, help="Accuracy of Douglas-Peucker line simplification, measured in degrees. Lower values give more nodes, and more accurate lines. Defaults to 0.0003.")
|
|---|
| 439 | parser.add_option("--dp-nr", action="store_true", dest="dp_nr", default=False, help="Use experimental non-recursive DP implementation")
|
|---|
| 440 | parser.add_option("--vr", action="store_true", dest="use_vr", default=False, help="Use vertex reduction before applying line simplification (off by default).")
|
|---|
| 441 | parser.add_option("--vr-epsilon", type="float", default=0.0005, metavar="RADIUS", help="Radius used for vertex reduction (measured in degrees). Defaults to 0.0005.")
|
|---|
| 442 | parser.add_option("--water", action="store_const", const="water", dest="natural_type", default="coastline", help="Tag ways as natural=water instead of natural=coastline")
|
|---|
| 443 | parser.add_option("--left", type="float", metavar="LONGITUDE", default=-180, help="Left (west) longitude for bounding box")
|
|---|
| 444 | parser.add_option("--right", type="float", metavar="LONGITUDE", default=180, help="Right (east) longitude for bounding box")
|
|---|
| 445 | parser.add_option("--top", type="float", metavar="LATITUDE", default=90, help="Top (north) latitude for bounding box")
|
|---|
| 446 | parser.add_option("--bottom", type="float", metavar="LATITUDE", default=-90, help="Bottom (south) latitude for bounding box")
|
|---|
| 447 | parser.add_option("--josm", action="store_true", dest="josm_mode", default=False, help="Operate in JOSM plugin mode")
|
|---|
| 448 | parser.add_option("--wms", type="string", dest="wmslayer", default="IR1", help="WMS layer to use")
|
|---|
| 449 |
|
|---|
| 450 | global options # Ugly, I know...
|
|---|
| 451 | (options, args) = parser.parse_args()
|
|---|
| 452 |
|
|---|
| 453 | if len(args) > 0:
|
|---|
| 454 | parser.print_help()
|
|---|
| 455 | return
|
|---|
| 456 |
|
|---|
| 457 | (start_lat, start_lon, infile) = (options.lat, options.lon, options.infile)
|
|---|
| 458 |
|
|---|
| 459 | if (start_lat is None or start_lon is None) and infile is None:
|
|---|
| 460 | if not options.josm_mode:
|
|---|
| 461 | parser.print_help()
|
|---|
| 462 | print
|
|---|
| 463 | error("Error: you must specify a starting latitude and longitude.")
|
|---|
| 464 | return
|
|---|
| 465 |
|
|---|
| 466 | if infile is not None:
|
|---|
| 467 | if start_lat is not None or start_lon is not None:
|
|---|
| 468 | error("Error: you cannot use both --infile and --lat or --lon.")
|
|---|
| 469 | return
|
|---|
| 470 | try:
|
|---|
| 471 | locs = get_locs(infile)
|
|---|
| 472 | #print locs
|
|---|
| 473 | except FatalError, e:
|
|---|
| 474 | error("%s" % e)
|
|---|
| 475 | return
|
|---|
| 476 | else:
|
|---|
| 477 | locs = [((start_lat, start_lon), options.threshold)]
|
|---|
| 478 |
|
|---|
| 479 | dirname = options.startdir.lower()
|
|---|
| 480 | if dirname in dirnames:
|
|---|
| 481 | startdir = dirnames.index(dirname)
|
|---|
| 482 | elif dirname in dirabbrevs:
|
|---|
| 483 | startdir = dirabbrevs.index(dirname)
|
|---|
| 484 | else:
|
|---|
| 485 | error("Error: Can't understand starting direction \"%s\". Vaild options are %s." % (dirname, ", ".join(dirnames + dirabbrevs)))
|
|---|
| 486 | return
|
|---|
| 487 |
|
|---|
| 488 | message("Starting direction is %s." % dirnames[startdir])
|
|---|
| 489 |
|
|---|
| 490 | bbox = BBox(options.top, options.left, options.bottom, options.right)
|
|---|
| 491 |
|
|---|
| 492 | try:
|
|---|
| 493 | lakes = []
|
|---|
| 494 | for (loc, threshold) in locs:
|
|---|
| 495 | nodes = trace_lake(loc, threshold, startdir, bbox)
|
|---|
| 496 |
|
|---|
| 497 | message("%d nodes generated." % len(nodes))
|
|---|
| 498 |
|
|---|
| 499 | if len(nodes) > 0:
|
|---|
| 500 | if options.use_vr:
|
|---|
| 501 | nodes = vertex_reduce(nodes, options.vr_epsilon)
|
|---|
| 502 | message("After vertex reduction, %d nodes remain." % len(nodes))
|
|---|
| 503 |
|
|---|
| 504 | if options.use_dp:
|
|---|
| 505 | try:
|
|---|
| 506 | if options.dp_nr:
|
|---|
| 507 | nodes = douglas_peucker_nonrecursive(nodes, options.dp_epsilon)
|
|---|
| 508 | #print "Final result: %s" % (nodes,)
|
|---|
| 509 | else:
|
|---|
| 510 | nodes = douglas_peucker(nodes, options.dp_epsilon)
|
|---|
| 511 | message("After Douglas-Peucker approximation, %d nodes remain." % len(nodes))
|
|---|
| 512 | except FatalError, e:
|
|---|
| 513 | raise e
|
|---|
| 514 | except:
|
|---|
| 515 | raise FatalError("Line simplification failed -- there are probably too many nodes.")
|
|---|
| 516 |
|
|---|
| 517 | lakes.append(nodes)
|
|---|
| 518 |
|
|---|
| 519 | if options.josm_mode:
|
|---|
| 520 | output_to_josm(lakes)
|
|---|
| 521 | else:
|
|---|
| 522 | message("Writing to %s" % options.outfile)
|
|---|
| 523 | f = open(options.outfile, "w")
|
|---|
| 524 | write_osm(f, lakes, options.waylength)
|
|---|
| 525 | f.close()
|
|---|
| 526 |
|
|---|
| 527 | #tiles = tilelist(nodes)
|
|---|
| 528 |
|
|---|
| 529 | #for tile in tiles:
|
|---|
| 530 | # print "./tilesGen.pl xy %d %d; ./upload.pl" % tile
|
|---|
| 531 |
|
|---|
| 532 | #print tiles
|
|---|
| 533 | except FatalError, e:
|
|---|
| 534 | error("%s" % e)
|
|---|
| 535 | error("Bailing out...")
|
|---|
| 536 |
|
|---|
| 537 | if __name__ == "__main__":
|
|---|
| 538 | main()
|
|---|
| 539 |
|
|---|
