OpenStreetMap logo OpenStreetMap

h4ck3rm1k3's Comments

Post When Comment
seextractor on osm

LOL
seextractor on osm

well, of course you need to swap lat and long... doh.
seextractor on osm

## copyright 2009 james michael dupont
## AGPL 3.0
## converts the CATS results into OSM format.
## example usage :
# convert landsat_4000_2000_xy_86000_168000.png -depth 16 landsat_4000_2000_xy_86000_168000.fits
# sextractor landsat_4000_2000_xy_86000_168000.fits
# perl ./convert.pl 86000 168000 test.cat > test.osm

use strict;
use warnings;

# bottom_left
my $basex = shift @ARGV;
my $basey = shift @ARGV;
my $filename = shift @ARGV;

my $tilesize = 2000;

my $topleft = $basey - $tilesize;

my $resolution = 4000;

print '';

my $id = 0;
open IN, $filename;
while ()
{
next if /\#/;
my ($test,$X_IMAGE,$Y_IMAGE,$X2,$Y2) = split (/\s+/);

print "\n";

$id --;
# print join ("!",($NUMBER,$FLUXERR_ISO,$FLUX_AUTO,$FLUXERR_AUTO,$X_IMAGE,$Y_IMAGE,$FLAGS)) . "\n";

my $newy = ($topleft + $Y_IMAGE) / $resolution;
my $newx = ($basex + $X_IMAGE) / $resolution;
my $newy2 = ($topleft + $Y2) / $resolution;
my $newx2 = ($basex + $X2) / $resolution;

# print "Data:" . $Y_IMAGE . "/" . $X_IMAGE . "\n";
print "\n";
$id --;
print "\n";
$id --;
print "\n";
$id --;
print "\n";

$id --;

my $id1 = $id + 4;
my $id2 = $id + 3;
my $id3 = $id + 2;
my $id4 = $id + 1;
print "








";

#print $newy . "/" . $newx . "\n";

}

print "";
close IN;
remove duplicate points from osm

1. I need to understand the problem first.
2. Forestwalker takes a while to run for such a huge input, and if it works, why should I break it.
when i have solved the problem, I will look at fixing forest walker.

merge multiple files mergemultiple.pl

Thanks, I will look at that. I know osmosis is powerful.
ForestWalker.py

prizren> i use the diary as my source control
let me walk you through it
first you need these files
http://svn.openstreetmap.org/applications/utils/import/lakewalker/OSMReader.py
http://svn.openstreetmap.org/applications/utils/import/lakewalker/OSM.py
then you need this file
http://www.pastebin.ca/1532240
and then it takes parameters of a osm file
or just a point
and the threshold
python ./lakewalker.py --lat 42.4052 --lon 20.4782 --threshold 80 -o forest.osm
or
python ./lakewalker.py -i somedots.osm --threshold 80 -o forest.osm

ForestWalker.py

here is my new hack that will create a grid of all the points in a 

#!/usr/bin/python



# Forestwalker 0.0003-gridhack.

# 2007-07-09: Initial public release (v0.2) by Dshpak

# 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().

# 2007-08-04: Added experimental non-recursive Douglas-Peucker algorithm (--dp-nr). Added --startdir option.

# 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.

# 2009-08-17 James Michael DuPont h4ck3rm1k3@flossk.org working



"""Lakewalker II - A tool to automatically download and trace Landsat imagery, to generate OpenStreetMap data.



Requires the Python Imaging Library, available from http://www.pythonware.com/products/pil/"""



version="Lakewalker II v0.6"



# TODO:

# - Accept threshold tags in OSM input files.

# - Command line options:

#   - direction (deosil/widdershins)

#   - layer to use (monochrome only?)

#   - additional tags for the way

#   - Landsat download retries (count and delay)

#   - radii for loop detection

#   - fname for z12 tile output (list or script)

#   - path to tilesGen for z12 script output

#   - debug mode (extra tags on nodes) (make this non-uploadable?)

# - The "got stuck" message should output coords

# - Better "got stuck" detection (detect mini loops)

# - Offset nodes outwards by a half-pixel or so, to prevent duplicate segments

# - Automatic threshold detection

# - (Correct/tested) non-recursive DP simplification

#

#

# For JOSM integration:

# - Add a --tmpdir option that controls the location of the Landsat

#   tiles, the results text file, and the lake.osm file (unless it's

#   got an absolute path

# - Add a --nocache option that keeps WMS tiles in memory but not on disk.



import math

import os

import urllib

from PIL import Image

import OSMReader

import time

import optparse



options = None



start_radius_big = 0.001

start_radius_small = 0.0002



## for tile walking

stepsize = 33



#dirs = ((1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1))

# steps of 10.. walk faster



stepsize = 1

errortries=20



dirs = ((1 * stepsize, 0), (1 * stepsize, 1 * stepsize), (0, 1 * stepsize), (-1 * stepsize, 1 * stepsize), (-1 * stepsize, 0), (-1 * stepsize, -1 * stepsize), (0, -1 * stepsize), (1 * stepsize, -1 * stepsize))



dirnames = ["east", "northeast", "north", "northwest", "west", "southwest", "south", "southeast"]

dirabbrevs = ["e", "ne", "n", "nw", "w", "sw", "s", "se"]



class FatalError(Exception):

    pass



def message(s):

    if options.josm_mode:

        print "m %s" % s

    else:

        print s

        

def error(s):

    if options.josm_mode:

        print "e %s" % s

    else:

        print s

        

class BBox:

    def __init__(self, top = 90, left = -180, bottom = -90, right = 180):

        self.left = left

        self.right = right

        self.top = top

        self.bottom = bottom

    def contains(self, loc):

        (lat, lon) = loc

        if lat > self.top or lat < self.bottom:

            return False

        if (self.right - self.left) % 360 == 0:

            return True

        return (lon - self.left) % 360 <= (self.right - self.left) % 360



def download_landsat(c1, c2, width, height, fname):

    layer = "global_mosaic_base"

#    style = "IR1"

    style = "IR2"



    (min_lat, min_lon) = c1

    (max_lat, max_lon) = c2

    

    message("Downloading Landsat tile for (%.6f,%.6f)-(%.6f,%.6f)." % (min_lat, min_lon, max_lat, max_lon))

    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)

    print url

    try:

        urllib.urlretrieve(url, fname)

    except IOError, e:

        raise FatalError("Error downloading tile: %s" % e.strerror)

    if not os.path.exists(fname):

        raise FatalError("Error: Could not retreive url %s" % url)

    print "Got it %s" % fname



def xy_to_geo(xy):

    (x, y) = xy

    (lat, lon) = (y / float(options.resolution), x / float(options.resolution))

    return (lat, lon)



def geo_to_xy(geo):

    (lat, lon) = geo

    coord = lambda L: math.floor(L * options.resolution + 0.5)

    (x, y) = (coord(lon), coord(lat))

    return (x, y)



class WMSManager:

    def __init__(self):

        self.images = {}



    def get_tile(self, xy):



        im = None

        while im is None :

            (x, y) = xy

            bottom_left_xy = (int(math.floor(x / options.tilesize)) * options.tilesize, int(math.floor(y / options.tilesize)) * options.tilesize)

            top_right_xy = (bottom_left_xy[0] + options.tilesize, bottom_left_xy[1] + options.tilesize)

            fname = "landsat_%d_%d_xy_%d_%d.png" % (options.resolution, options.tilesize, bottom_left_xy[0], bottom_left_xy[1])

            im = self.images.get(fname, None)



            if not os.path.exists(fname):

                raise FatalError("Error: Could not get image file %s" % fname)

            try:

                im = Image.open(fname)

                self.images[fname] = im

                message("Using imagery in %s for %s." % (fname, xy_to_geo(xy)))

            except IOError:

                error("Download was corrupt...Deleting %s..." % fname)

                raise FatalError("Error: Could not get image file %s" % fname)



        return (im, bottom_left_xy)

    

    def get_pixel(self, xy):

        (x, y) = xy

        (tile, (tx, ty)) = self.get_tile(xy)

        tile_xy = (x - tx, (options.tilesize - 1) - (y - ty))

        print "%s maps to %s" % (xy, tile_xy)

        return tile.getpixel(tile_xy)



def walk_tile(loc, threshold, start_dir, bbox):

    wms = WMSManager()

    xy = geo_to_xy(loc)

    nodelist = []

    message("Starting coordinate: %.4f, %.4f" % loc)

    message("Starting position: %d, %d" % xy)



    #v = wms.get_pixel(xy)

    #get the tile that the pixel contains

    (tile, (tx, ty)) = wms.get_tile(xy)

    



#    print tile

    for x in xrange(tx +1, tx + options.tilesize -1, stepsize) :

        for y in xrange(ty+1, ty + options.tilesize -1, stepsize) :

            loc = xy_to_geo((x, y))

            nodelist.append(loc)



#    message("Tile : %d, %d" % (tx, ty))

    return nodelist



def trace_lake(loc, threshold, start_dir, bbox):

    wms = WMSManager()

    xy = geo_to_xy(loc)

    nodelist = []



    message("Starting coordinate: %.4f, %.4f" % loc)

    message("Starting position: %d, %d" % xy)



    #v = wms.get_pixel(xy)

    #get the tile that the pixel contains

    (tile, (tx, ty)) = wms.get_tile(xy)



    

    if not bbox.contains(loc):

        raise FatalError("Error: Starting location is outside bounding box!")



    errcount = errortries



    while True:

        loc = xy_to_geo(xy)

        if not bbox.contains(loc):

            break

        

        v = wms.get_pixel(xy)

        if v > threshold:

            errcount = errcount -1 # allow for three errors

            if errcount < 0 :

                break



        xy = (xy[0] + dirs[start_dir][0], xy[1] + dirs[start_dir][1])



    start_xy = xy

    start_loc = xy_to_geo(xy)

    print xy 

    message("Found shore at  lat %.4f lon %.4f" % start_loc)



    #dirs = ((1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1), (0, 1), (1, 1))

    last_dir = start_dir



    detect_loop = False

    

    for i in xrange(options.maxnodes):

        if i % 250 == 0:

            if i > 0:

                message("%s nodes so far..." % i)



        for d in xrange(1, len(dirs)):

            new_dir = dirs[(last_dir + d + 4) % 8]

            test_xy = (xy[0] + new_dir[0], xy[1] + new_dir[1])

            test_loc = xy_to_geo(test_xy)

            if not bbox.contains(test_loc):

                break

            

            v = wms.get_pixel(test_xy)

            print "%s: %s: %s" % (new_dir, test_xy, v)

            if v > threshold:

                break



        if d == 8:

            error("Got stuck.")

            break



        print "Moving to %s, direction %s (was %s)" % (test_xy, new_dir, dirs[last_dir])

        last_dir = (last_dir + d + 4) % 8

        xy = test_xy



        if xy == start_xy:

            break



        loc = xy_to_geo(xy)

        nodelist.append(loc)

        

        start_proximity = (loc[0] - start_loc[0]) ** 2 + (loc[1] - start_loc[1]) ** 2

        if detect_loop:

            if start_proximity < start_radius_small ** 2:

                break

        else:

            if start_proximity > start_radius_big ** 2:

                detect_loop = True

    return nodelist



def vertex_reduce(nodes, proximity):

    test_v = nodes[0]

    reduced_nodes = [test_v]

    prox_sq = proximity ** 2

    for v in nodes:

        if (v[0] - test_v[0]) ** 2 + (v[1] - test_v[1]) ** 2 > prox_sq:

            reduced_nodes.append(v)

            test_v = v

    return reduced_nodes



def point_line_distance(p0, p1, p2):

    ((x0, y0), (x1, y1), (x2, y2)) = (p0, p1, p2)



    if x2 == x1 and y2 == y1:

        # Degenerate cast: the "line" is actually a point.

        return math.sqrt((x1-x0)**2 + (y1-y0)**2)

    else:

        # I don't understand this at all. Thank you, Mathworld.

        # http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html

        return abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / math.sqrt((x2-x1)**2 + (y2-y1)**2)



def douglas_peucker(nodes, epsilon):

    print "Running DP on %d nodes" % len(nodes)

    farthest_node = None

    farthest_dist = 0

    first = nodes[0]

    last = nodes[-1]



    for i in xrange(1, len(nodes) - 1):

        d = point_line_distance(nodes[i], first, last)

        if d > farthest_dist:

            farthest_dist = d

            farthest_node = i



    if farthest_dist > epsilon:

        seg_a = douglas_peucker(nodes[0:farthest_node+1], epsilon)

        seg_b = douglas_peucker(nodes[farthest_node:-1], epsilon)

        print "Minimized %d nodes to %d + %d nodes" % (len(nodes), len(seg_a), len(seg_b))

        nodes = seg_a[:-1] + seg_b

    else:

        return [nodes[0], nodes[-1]]



    return nodes



def dp_findpoint(nodes, start, end):

    farthest_node = None

    farthest_dist = 0

    print "dp_findpoint(nodes, %s, %s)" % (start, end)

    first = nodes[start]

    last = nodes[end]



    for i in xrange(start + 1, end):

        d = point_line_distance(nodes[i], first, last)

        if d > farthest_dist:

            farthest_dist = d

            farthest_node = i



    return (farthest_node, farthest_dist)



def douglas_peucker_nonrecursive(nodes, epsilon):

    print "Running DP on %d nodes" % len(nodes)

    command_stack = [(0, len(nodes) - 1)]

    result_stack = []



    while len(command_stack) > 0:

        cmd = command_stack.pop()

        if type(cmd) == tuple:

            (start, end) = cmd

            (node, dist) = dp_findpoint(nodes, start, end)

            if dist > epsilon:

                command_stack.append("+")

                command_stack.append((start, node))

                command_stack.append((node, end))

            else:

                result_stack.append((start, end))

        elif cmd == "+":

            first = result_stack.pop()

            second = result_stack.pop()

            if first[-1] == second[0]:

                result_stack.append(first + second[1:])

                print "Added %s and %s; result is %s" % (first, second, result_stack[-1])

            else:

                error("ERROR: Cannot connect nodestrings!")

                #print first

                #print second

                return

        else:

            error("ERROR: Can't understand command \"%s\"" % (cmd,))

            return



    if len(result_stack) == 1:

        return [nodes[x] for x in result_stack[0]]

    else:

        error("ERROR: Command stack is empty but result stack has %d nodes!" % len(result_stack))

        return



    farthest_node = None

    farthest_dist = 0

    first = nodes[0]

    last = nodes[-1]



    for i in xrange(1, len(nodes) - 1):

        d = point_line_distance(nodes[i], first, last)

        if d > farthest_dist:

            farthest_dist = d

            farthest_node = i



    if farthest_dist > epsilon:

        seg_a = douglas_peucker(nodes[0:farthest_node+1], epsilon)

        seg_b = douglas_peucker(nodes[farthest_node:-1], epsilon)

        print "Minimized %d nodes to %d + %d nodes" % (len(nodes), len(seg_a), len(seg_b))

        nodes = seg_a[:-1] + seg_b

    else:

        return [nodes[0], nodes[-1]]



    return nodes



def output_to_josm(lakelist):

    # Description of JOSM output format:

    # m text - Status message text, to be displayed in a display window

    # e text - Error message text

    # s nnn - Start full node list, nnn tracings following

    # 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)

    # n lat lon - A node

    # x - End of data

    print "s %s" % len(lakelist)

    for nodelist in lakelist:

        print "t %s" % len(nodelist)

        for node in nodelist:

            print "n %.7f %.7f" % (node[0], node[1])

    print "x"



def write_osm(f, lakelist, waysize):

    f.write('')

    cur_id = -1

    way_count = 0

    for nodelist in lakelist:

        first_node_id = cur_id

        cur_way = []

        ways = [cur_way]

        for loc in nodelist:

            #f.write('  \n' % (cur_id, loc[0], loc[1], -cur_id, geo_to_xy(loc)[0], geo_to_xy(loc)[1]))

            f.write('' % (cur_id, loc[0], loc[1]))

            last_node_id = cur_id

            cur_id = cur_id - 1



        print "Nodes: %d, %d" % (first_node_id, last_node_id)



#        f.write('\n' %  cur_id)



#        first_segment_id = cur_id



#        for seg in xrange(first_node_id, last_node_id, -1):

#            f.write('' % ( seg))

            #

#            cur_id = cur_id - 1

#        f.write('' % (cur_id, last_node_id, first_node_id))



#        f.write('\n')



#        last_segment_id = cur_id

        cur_id = cur_id - 1



        print "Segments: %d, %d" % (first_node_id, last_node_id)

        

        for seg in xrange(first_node_id, last_node_id - 1, -1):

            if len(cur_way) >= waysize:

                cur_way = []

                ways.append(cur_way)

            cur_way.append(seg)

        for way in ways:

            f.write('\n' % cur_id)

            cur_id = cur_id - 1

            for seg in way:

                f.write('' % seg)



            f.write('' % first_node_id) # close it

            f.write('' % options.natural_type)

            f.write('')

            f.write('\n')

        way_count = way_count + len(ways)

    

    f.write('')

    message("Generated %d %s." % (way_count, ["way", "ways"][way_count > 0]))



def outputlakes(lakes):



    filenum =0



    if options.josm_mode:

        output_to_josm(lakes)

    else:

        

        filename = "%s-%d.osm" % (options.outfile, filenum)        

        filenum = filenum +1 

        message("Writing to %s" % filename)

        f = open(filename, "w")

        write_osm(f, lakes, options.waylength)

        f.close()

        message("wrote the file %s" % filename)

        

        lakes = []

    return lakes



def get_locs(infile):

    nodes = []

    segments = []

    ways = []

    reader = OSMReader.OSMReader()

    reader.nodeHandler = lambda x: nodes.append(x)

    reader.segmentHandler = lambda x: segments.append(x)

    reader.wayHandler = lambda x: ways.append(x)

    reader.run(file(infile))

    if len(segments) > 0 or len(ways) > 0:

        raise FatalError("Error: Input file must only contain nodes -- no segments or ways.")

    if len(nodes) == 0:

        raise FatalError("Error: No nodes found in input file.")

    return [((node.lat, node.lon), options.threshold) for node in nodes]



def main():

    parser = optparse.OptionParser(version=version)



    parser.add_option("--lat", type="float", metavar="LATITUDE", help="Starting latitude. Required, unless --infile is used..")

    parser.add_option("--lon", type="float", metavar="LONGITUDE", help="Starting longitude. Required, unless --infile is used.")

    parser.add_option("--startdir", type="string", default="east", metavar="DIR", help="Direction to travel from start position when seeking land. Defaults to \"east\".")

    parser.add_option("--infile", "-i", type="string", metavar="FILE", help="OSM file containing nodes representing starting points.")

    parser.add_option("--out", "-o", default="forest.osm", dest="outfile", metavar="FILE", help="Output filename. Defaults to forest.osm.")

    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.")

    parser.add_option("--maxnodes", type="int", default="50000", metavar="N", help="Maximum number of nodes to generate before bailing out. Defaults to 50000.")

    parser.add_option("--waylength", type="int", default=250, metavar="MAXLEN", help="Maximum nuber of nodes allowed in one way. Defaults to 250.")

    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.")

    parser.add_option("--tilesize", type="int", default=2000, help="Size of one landsat tile, measured in pixels. Defaults to 2000.")

    parser.add_option("--no-dp", action="store_false", dest="use_dp", default=True, help="Disable Douglas-Peucker line simplification (not recommended)")

    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.")

    parser.add_option("--dp-nr", action="store_true", dest="dp_nr", default=False, help="Use experimental non-recursive DP implementation")

    parser.add_option("--vr", action="store_true", dest="use_vr", default=False, help="Use vertex reduction before applying line simplification (off by default).")

    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.")

#    parser.add_option("--water", action="store_const", const="water", dest="natural_type", default="coastline", help="Tag ways as natural=water instead of natural=coastline")

    parser.add_option("--forest", action="store_const", const="woords", dest="natural_type", default="forest", help="Tag ways as natural=woods instead of natural=forest")

    #options.natural_type



    parser.add_option("--left", type="float", metavar="LONGITUDE", default=-180, help="Left (west) longitude for bounding box")

    parser.add_option("--right", type="float", metavar="LONGITUDE", default=180, help="Right (east) longitude for bounding box")

    parser.add_option("--top", type="float", metavar="LATITUDE", default=90, help="Top (north) latitude for bounding box")

    parser.add_option("--bottom", type="float", metavar="LATITUDE", default=-90, help="Bottom (south) latitude for bounding box")

    parser.add_option("--josm", action="store_true", dest="josm_mode", default=False, help="Operate in JOSM plugin mode")



    global options # Ugly, I know...

    (options, args) = parser.parse_args()



    if len(args) > 0:

        parser.print_help()

        return



    (start_lat, start_lon, infile) = (options.lat, options.lon, options.infile)



    if (start_lat is None or start_lon is None) and infile is None:

        if not options.josm_mode:

            parser.print_help()

            print

        error("Error: you must specify a starting latitude and longitude.")

        return



    if infile is not None:

        if start_lat is not None or start_lon is not None:

            error("Error: you cannot use both --infile and --lat or --lon.")

            return

        try:

            locs = get_locs(infile)

            #print locs

        except FatalError, e:

            error("%s" % e)

            return

    else:

        locs = [((start_lat, start_lon), options.threshold)]



    dirname = options.startdir.lower()

    if dirname in dirnames:

        startdir = dirnames.index(dirname)

    elif dirname in dirabbrevs:

        startdir = dirabbrevs.index(dirname)

    else:

        error("Error: Can't understand starting direction \"%s\". Vaild options are %s." % (dirname, ", ".join(dirnames + dirabbrevs)))

        return



    message("Starting direction is %s." % dirnames[startdir])



    bbox = BBox(options.top, options.left, options.bottom, options.right)



    try:

        lakes = []

        filenum=1



        # create a grid for just the first point

        print locs[0][0]

        locgrid = walk_tile(locs[0][0], 0, startdir, bbox)



        nodes = []



        for (loc) in locgrid:

            #threshold

            #(lat, lon) = loc

            #newloc = (x, y) = (coord(lon), coord(lat))            



            try:

                print "tracing at " 

                print loc



                #nodes2 = trace_lake(loc, options.threshold, startdir, bbox)

                nodes.append(loc)

            except FatalError, e:

                error("%s" % e)



            message("%d nodes generated." % len(nodes))



            # skip over very small groups

            if len(nodes) > 0:

                if options.use_vr:

                    nodes = vertex_reduce(nodes, options.vr_epsilon)

                    message("After vertex reduction, %d nodes remain." % len(nodes))



                if options.use_dp:

                    try:

                        if options.dp_nr:

                            nodes = douglas_peucker_nonrecursive(nodes, options.dp_epsilon)

                            #print "Final result: %s" % (nodes,)

                        else:

                            nodes = douglas_peucker(nodes, options.dp_epsilon)

                        message("After Douglas-Peucker approximation, %d nodes remain." % len(nodes))

                    except FatalError, e:

                        raise e

                    except:

                        message(" Douglas-Peucker approximation failed, who cares., %d nodes remain." % len(nodes))

                        #raise FatalError("Line simplification failed -- there are probably too many nodes.")



                lakes.append(nodes)



            lakes= outputlakes(lakes)



    except FatalError, e:

        error("%s" % e)

        error("Bailing out...")



if __name__ == "__main__":

    main()

ForestWalker.py

Sorry, you need the other files here :
http://svn.openstreetmap.org/applications/utils/import/lakewalker/
http://svn.openstreetmap.org/applications/utils/import/lakewalker/OSMReader.py
http://svn.openstreetmap.org/applications/utils/import/lakewalker/OSM.py
ForestWalker.py

Original source
http://svn.openstreetmap.org/applications/utils/import/lakewalker/lakewalker.py

JOSM and Pen Tablet

I got it
osm.org/user/h4ck3rm1k3/diary/7538
Forest walker.
ForestWalker.py 

#!/usr/bin/python



# Forestwalker 0.0001.

# 2007-07-09: Initial public release (v0.2) by Dshpak

# 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().

# 2007-08-04: Added experimental non-recursive Douglas-Peucker algorithm (--dp-nr). Added --startdir option.

# 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.

# 2009-08-17 James Michael DuPont h4ck3rm1k3@flossk.org working



"""Lakewalker II - A tool to automatically download and trace Landsat imagery, to generate OpenStreetMap data.



Requires the Python Imaging Library, available from http://www.pythonware.com/products/pil/"""



version="Lakewalker II v0.6"



# TODO:

# - Accept threshold tags in OSM input files.

# - Command line options:

#   - direction (deosil/widdershins)

#   - layer to use (monochrome only?)

#   - additional tags for the way

#   - Landsat download retries (count and delay)

#   - radii for loop detection

#   - fname for z12 tile output (list or script)

#   - path to tilesGen for z12 script output

#   - debug mode (extra tags on nodes) (make this non-uploadable?)

# - The "got stuck" message should output coords

# - Better "got stuck" detection (detect mini loops)

# - Offset nodes outwards by a half-pixel or so, to prevent duplicate segments

# - Automatic threshold detection

# - (Correct/tested) non-recursive DP simplification

#

#

# For JOSM integration:

# - Add a --tmpdir option that controls the location of the Landsat

#   tiles, the results text file, and the lake.osm file (unless it's

#   got an absolute path

# - Add a --nocache option that keeps WMS tiles in memory but not on disk.



import math

import os

import urllib

from PIL import Image

import OSMReader

import time

import optparse



options = None



start_radius_big = 0.001

start_radius_small = 0.0002



dirs = ((1, 0), (1, 1), (0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1))

dirnames = ["east", "northeast", "north", "northwest", "west", "southwest", "south", "southeast"]

dirabbrevs = ["e", "ne", "n", "nw", "w", "sw", "s", "se"]



class FatalError(Exception):

    pass



def message(s):

    if options.josm_mode:

        print "m %s" % s

    else:

        print s

        

def error(s):

    if options.josm_mode:

        print "e %s" % s

    else:

        print s

        

class BBox:

    def __init__(self, top = 90, left = -180, bottom = -90, right = 180):

        self.left = left

        self.right = right

        self.top = top

        self.bottom = bottom

    def contains(self, loc):

        (lat, lon) = loc

        if lat > self.top or lat < self.bottom:

            return False

        if (self.right - self.left) % 360 == 0:

            return True

        return (lon - self.left) % 360 <= (self.right - self.left) % 360



def download_landsat(c1, c2, width, height, fname):

    layer = "global_mosaic_base"

#    style = "IR1"

    style = "IR2"



    (min_lat, min_lon) = c1

    (max_lat, max_lon) = c2

    

    message("Downloading Landsat tile for (%.6f,%.6f)-(%.6f,%.6f)." % (min_lat, min_lon, max_lat, max_lon))

    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)

    print url

    try:

        urllib.urlretrieve(url, fname)

    except IOError, e:

        raise FatalError("Error downloading tile: %s" % e.strerror)

    if not os.path.exists(fname):

        raise FatalError("Error: Could not retreive url %s" % url)

    print "Got it %s" % fname



def xy_to_geo(xy):

    (x, y) = xy

    (lat, lon) = (y / float(options.resolution), x / float(options.resolution))

    return (lat, lon)



def geo_to_xy(geo):

    (lat, lon) = geo

    coord = lambda L: math.floor(L * options.resolution + 0.5)

    (x, y) = (coord(lon), coord(lat))

    return (x, y)



class WMSManager:

    def __init__(self):

        self.images = {}



    def get_tile(self, xy):

        fail_count = 0

        im = None

        while im is None and fail_count < 4:

            (x, y) = xy

            bottom_left_xy = (int(math.floor(x / options.tilesize)) * options.tilesize, int(math.floor(y / options.tilesize)) * options.tilesize)

            top_right_xy = (bottom_left_xy[0] + options.tilesize, bottom_left_xy[1] + options.tilesize)

            fname = "landsat_%d_%d_xy_%d_%d.png" % (options.resolution, options.tilesize, bottom_left_xy[0], bottom_left_xy[1])

            im = self.images.get(fname, None)

            if im is None:

                if not os.path.exists(fname):

                    bottom_left = xy_to_geo(bottom_left_xy)

                    top_right = xy_to_geo(top_right_xy)

                    download_landsat(bottom_left, top_right, options.tilesize, options.tilesize, fname)

                if not os.path.exists(fname):

                    raise FatalError("Error: Could not get image file %s" % fname)

                try:

                    im = Image.open(fname)

                    self.images[fname] = im

                    message("Using imagery in %s for %s." % (fname, xy_to_geo(xy)))

                except IOError:

                    error("Download was corrupt...Deleting %s..." % fname)

                    #os.unlink(fname)

                    im = None



                if im is None:

                    message("Sleeping and retrying download...")

                    time.sleep(4)

                    fail_count = fail_count + 1



        if im is None:

            #if os.path.exists(fname):

            #    print open(fname).readlines()

            raise FatalError("Couldn't get image file %s." % fname)



        #return (im, top_left_xy)

        return (im, bottom_left_xy)

    

    def get_pixel(self, xy):

        (x, y) = xy

        (tile, (tx, ty)) = self.get_tile(xy)

        tile_xy = (x - tx, (options.tilesize - 1) - (y - ty))

        print "%s maps to %s" % (xy, tile_xy)

        return tile.getpixel(tile_xy)



def trace_lake(loc, threshold, start_dir, bbox):

    wms = WMSManager()

    xy = geo_to_xy(loc)

    nodelist = []



    message("Starting coordinate: %.4f, %.4f" % loc)

    message("Starting position: %d, %d" % xy)



    if not bbox.contains(loc):

        raise FatalError("Error: Starting location is outside bounding box!")



    while True:

        loc = xy_to_geo(xy)

        if not bbox.contains(loc):

            break

        

        v = wms.get_pixel(xy)

        if v > threshold:

            break



        xy = (xy[0] + dirs[start_dir][0], xy[1] + dirs[start_dir][1])



    start_xy = xy

    start_loc = xy_to_geo(xy)

    print xy 

    message("Found shore at  lat %.4f lon %.4f" % start_loc)



    #dirs = ((1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1), (0, 1), (1, 1))

    last_dir = start_dir



    detect_loop = False

    

    for i in xrange(options.maxnodes):

        if i % 250 == 0:

            if i > 0:

                message("%s nodes so far..." % i)



        for d in xrange(1, len(dirs)):

            new_dir = dirs[(last_dir + d + 4) % 8]

            test_xy = (xy[0] + new_dir[0], xy[1] + new_dir[1])

            test_loc = xy_to_geo(test_xy)

            if not bbox.contains(test_loc):

                break

            

            v = wms.get_pixel(test_xy)

            print "%s: %s: %s" % (new_dir, test_xy, v)

            if v > threshold:

                break



        if d == 8:

            error("Got stuck.")

            break



        print "Moving to %s, direction %s (was %s)" % (test_xy, new_dir, dirs[last_dir])

        last_dir = (last_dir + d + 4) % 8

        xy = test_xy



        if xy == start_xy:

            break



        loc = xy_to_geo(xy)

        nodelist.append(loc)

        

        start_proximity = (loc[0] - start_loc[0]) ** 2 + (loc[1] - start_loc[1]) ** 2

        if detect_loop:

            if start_proximity < start_radius_small ** 2:

                break

        else:

            if start_proximity > start_radius_big ** 2:

                detect_loop = True

    return nodelist



def vertex_reduce(nodes, proximity):

    test_v = nodes[0]

    reduced_nodes = [test_v]

    prox_sq = proximity ** 2

    for v in nodes:

        if (v[0] - test_v[0]) ** 2 + (v[1] - test_v[1]) ** 2 > prox_sq:

            reduced_nodes.append(v)

            test_v = v

    return reduced_nodes



def point_line_distance(p0, p1, p2):

    ((x0, y0), (x1, y1), (x2, y2)) = (p0, p1, p2)



    if x2 == x1 and y2 == y1:

        # Degenerate cast: the "line" is actually a point.

        return math.sqrt((x1-x0)**2 + (y1-y0)**2)

    else:

        # I don't understand this at all. Thank you, Mathworld.

        # http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html

        return abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / math.sqrt((x2-x1)**2 + (y2-y1)**2)



def douglas_peucker(nodes, epsilon):

    print "Running DP on %d nodes" % len(nodes)

    farthest_node = None

    farthest_dist = 0

    first = nodes[0]

    last = nodes[-1]



    for i in xrange(1, len(nodes) - 1):

        d = point_line_distance(nodes[i], first, last)

        if d > farthest_dist:

            farthest_dist = d

            farthest_node = i



    if farthest_dist > epsilon:

        seg_a = douglas_peucker(nodes[0:farthest_node+1], epsilon)

        seg_b = douglas_peucker(nodes[farthest_node:-1], epsilon)

        print "Minimized %d nodes to %d + %d nodes" % (len(nodes), len(seg_a), len(seg_b))

        nodes = seg_a[:-1] + seg_b

    else:

        return [nodes[0], nodes[-1]]



    return nodes



def dp_findpoint(nodes, start, end):

    farthest_node = None

    farthest_dist = 0

    print "dp_findpoint(nodes, %s, %s)" % (start, end)

    first = nodes[start]

    last = nodes[end]



    for i in xrange(start + 1, end):

        d = point_line_distance(nodes[i], first, last)

        if d > farthest_dist:

            farthest_dist = d

            farthest_node = i



    return (farthest_node, farthest_dist)



def douglas_peucker_nonrecursive(nodes, epsilon):

    print "Running DP on %d nodes" % len(nodes)

    command_stack = [(0, len(nodes) - 1)]

    result_stack = []



    while len(command_stack) > 0:

        cmd = command_stack.pop()

        if type(cmd) == tuple:

            (start, end) = cmd

            (node, dist) = dp_findpoint(nodes, start, end)

            if dist > epsilon:

                command_stack.append("+")

                command_stack.append((start, node))

                command_stack.append((node, end))

            else:

                result_stack.append((start, end))

        elif cmd == "+":

            first = result_stack.pop()

            second = result_stack.pop()

            if first[-1] == second[0]:

                result_stack.append(first + second[1:])

                print "Added %s and %s; result is %s" % (first, second, result_stack[-1])

            else:

                error("ERROR: Cannot connect nodestrings!")

                #print first

                #print second

                return

        else:

            error("ERROR: Can't understand command \"%s\"" % (cmd,))

            return



    if len(result_stack) == 1:

        return [nodes[x] for x in result_stack[0]]

    else:

        error("ERROR: Command stack is empty but result stack has %d nodes!" % len(result_stack))

        return



    farthest_node = None

    farthest_dist = 0

    first = nodes[0]

    last = nodes[-1]



    for i in xrange(1, len(nodes) - 1):

        d = point_line_distance(nodes[i], first, last)

        if d > farthest_dist:

            farthest_dist = d

            farthest_node = i



    if farthest_dist > epsilon:

        seg_a = douglas_peucker(nodes[0:farthest_node+1], epsilon)

        seg_b = douglas_peucker(nodes[farthest_node:-1], epsilon)

        print "Minimized %d nodes to %d + %d nodes" % (len(nodes), len(seg_a), len(seg_b))

        nodes = seg_a[:-1] + seg_b

    else:

        return [nodes[0], nodes[-1]]



    return nodes



def output_to_josm(lakelist):

    # Description of JOSM output format:

    # m text - Status message text, to be displayed in a display window

    # e text - Error message text

    # s nnn - Start full node list, nnn tracings following

    # 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)

    # n lat lon - A node

    # x - End of data

    print "s %s" % len(lakelist)

    for nodelist in lakelist:

        print "t %s" % len(nodelist)

        for node in nodelist:

            print "n %.7f %.7f" % (node[0], node[1])

    print "x"



def write_osm(f, lakelist, waysize):

    f.write('')

    cur_id = -1

    way_count = 0

    for nodelist in lakelist:

        first_node_id = cur_id

        cur_way = []

        ways = [cur_way]

        for loc in nodelist:

            #f.write('  \n' % (cur_id, loc[0], loc[1], -cur_id, geo_to_xy(loc)[0], geo_to_xy(loc)[1]))

            f.write('' % (cur_id, loc[0], loc[1]))

            last_node_id = cur_id

            cur_id = cur_id - 1



        print "Nodes: %d, %d" % (first_node_id, last_node_id)



#        f.write('\n' %  cur_id)



#        first_segment_id = cur_id



#        for seg in xrange(first_node_id, last_node_id, -1):

#            f.write('' % ( seg))

            #

#            cur_id = cur_id - 1

#        f.write('' % (cur_id, last_node_id, first_node_id))



#        f.write('\n')



#        last_segment_id = cur_id

        cur_id = cur_id - 1



        print "Segments: %d, %d" % (first_node_id, last_node_id)

        

        for seg in xrange(first_node_id, last_node_id - 1, -1):

            if len(cur_way) >= waysize:

                cur_way = []

                ways.append(cur_way)

            cur_way.append(seg)

        for way in ways:

            f.write('\n' % cur_id)

            cur_id = cur_id - 1

            for seg in way:

                f.write('' % seg)



            f.write('' % first_node_id) # close it

            f.write('' % options.natural_type)

            f.write('')

            f.write('\n')

        way_count = way_count + len(ways)

    

    f.write('')

    message("Generated %d %s." % (way_count, ["way", "ways"][way_count > 0]))



def get_locs(infile):

    nodes = []

    segments = []

    ways = []

    reader = OSMReader.OSMReader()

    reader.nodeHandler = lambda x: nodes.append(x)

    reader.segmentHandler = lambda x: segments.append(x)

    reader.wayHandler = lambda x: ways.append(x)

    reader.run(file(infile))

    if len(segments) > 0 or len(ways) > 0:

        raise FatalError("Error: Input file must only contain nodes -- no segments or ways.")

    if len(nodes) == 0:

        raise FatalError("Error: No nodes found in input file.")

    return [((node.lat, node.lon), options.threshold) for node in nodes]



def main():

    parser = optparse.OptionParser(version=version)



    parser.add_option("--lat", type="float", metavar="LATITUDE", help="Starting latitude. Required, unless --infile is used..")

    parser.add_option("--lon", type="float", metavar="LONGITUDE", help="Starting longitude. Required, unless --infile is used.")

    parser.add_option("--startdir", type="string", default="east", metavar="DIR", help="Direction to travel from start position when seeking land. Defaults to \"east\".")

    parser.add_option("--infile", "-i", type="string", metavar="FILE", help="OSM file containing nodes representing starting points.")

    parser.add_option("--out", "-o", default="forest.osm", dest="outfile", metavar="FILE", help="Output filename. Defaults to forest.osm.")

    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.")

    parser.add_option("--maxnodes", type="int", default="50000", metavar="N", help="Maximum number of nodes to generate before bailing out. Defaults to 50000.")

    parser.add_option("--waylength", type="int", default=250, metavar="MAXLEN", help="Maximum nuber of nodes allowed in one way. Defaults to 250.")

    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.")

    parser.add_option("--tilesize", type="int", default=2000, help="Size of one landsat tile, measured in pixels. Defaults to 2000.")

    parser.add_option("--no-dp", action="store_false", dest="use_dp", default=True, help="Disable Douglas-Peucker line simplification (not recommended)")

    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.")

    parser.add_option("--dp-nr", action="store_true", dest="dp_nr", default=False, help="Use experimental non-recursive DP implementation")

    parser.add_option("--vr", action="store_true", dest="use_vr", default=False, help="Use vertex reduction before applying line simplification (off by default).")

    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.")

#    parser.add_option("--water", action="store_const", const="water", dest="natural_type", default="coastline", help="Tag ways as natural=water instead of natural=coastline")

    parser.add_option("--forest", action="store_const", const="woords", dest="natural_type", default="forest", help="Tag ways as natural=woods instead of natural=forest")

    #options.natural_type



    parser.add_option("--left", type="float", metavar="LONGITUDE", default=-180, help="Left (west) longitude for bounding box")

    parser.add_option("--right", type="float", metavar="LONGITUDE", default=180, help="Right (east) longitude for bounding box")

    parser.add_option("--top", type="float", metavar="LATITUDE", default=90, help="Top (north) latitude for bounding box")

    parser.add_option("--bottom", type="float", metavar="LATITUDE", default=-90, help="Bottom (south) latitude for bounding box")

    parser.add_option("--josm", action="store_true", dest="josm_mode", default=False, help="Operate in JOSM plugin mode")



    global options # Ugly, I know...

    (options, args) = parser.parse_args()



    if len(args) > 0:

        parser.print_help()

        return



    (start_lat, start_lon, infile) = (options.lat, options.lon, options.infile)



    if (start_lat is None or start_lon is None) and infile is None:

        if not options.josm_mode:

            parser.print_help()

            print

        error("Error: you must specify a starting latitude and longitude.")

        return



    if infile is not None:

        if start_lat is not None or start_lon is not None:

            error("Error: you cannot use both --infile and --lat or --lon.")

            return

        try:

            locs = get_locs(infile)

            #print locs

        except FatalError, e:

            error("%s" % e)

            return

    else:

        locs = [((start_lat, start_lon), options.threshold)]



    dirname = options.startdir.lower()

    if dirname in dirnames:

        startdir = dirnames.index(dirname)

    elif dirname in dirabbrevs:

        startdir = dirabbrevs.index(dirname)

    else:

        error("Error: Can't understand starting direction \"%s\". Vaild options are %s." % (dirname, ", ".join(dirnames + dirabbrevs)))

        return



    message("Starting direction is %s." % dirnames[startdir])



    bbox = BBox(options.top, options.left, options.bottom, options.right)



    try:

        lakes = []

        for (loc, threshold) in locs:

            nodes = trace_lake(loc, threshold, startdir, bbox)



            message("%d nodes generated." % len(nodes))



            if len(nodes) > 0:

                if options.use_vr:

                    nodes = vertex_reduce(nodes, options.vr_epsilon)

                    message("After vertex reduction, %d nodes remain." % len(nodes))



                if options.use_dp:

                    try:

                        if options.dp_nr:

                            nodes = douglas_peucker_nonrecursive(nodes, options.dp_epsilon)

                            #print "Final result: %s" % (nodes,)

                        else:

                            nodes = douglas_peucker(nodes, options.dp_epsilon)

                        message("After Douglas-Peucker approximation, %d nodes remain." % len(nodes))

                    except FatalError, e:

                        raise e

                    except:

                        raise FatalError("Line simplification failed -- there are probably too many nodes.")



                lakes.append(nodes)



        if options.josm_mode:

            output_to_josm(lakes)

        else:

            message("Writing to %s" % options.outfile)

            f = open(options.outfile, "w")

            write_osm(f, lakes, options.waylength)

            f.close()

        

        #tiles = tilelist(nodes)

        

        #for tile in tiles:

        #    print "./tilesGen.pl xy %d %d; ./upload.pl" % tile

        

        #print tiles

    except FatalError, e:

        error("%s" % e)

        error("Bailing out...")



if __name__ == "__main__":

    main()

Mappping Party Prizren

http://www.google.com/calendar/event?action=VIEW&eid=bGh1azl0a2wwMnBmYjZyN3MxMWs2dWJrOTQgZnJlZS1zb2Z0d2FyZS1jb25mZXJlbmNlQGdvb2dsZWdyb3Vwcy5jb20&tok=MzAjamFtZXNtaWtlZHVwb250QGdvb2dsZW1haWwuY29tN2YwMWM1NmQwMzY5MDllYjJlMDRmNjFiMWY2NDI3NGUwNjVkYWEwZg&ctz=Europe%2FBerlin&hl=en&gsessionid=pvUok869GImuyCx92LR0Qw
Archeological Site map for Prizren/ Kosovo

Point 1 should be this church
osm.org/browse/node/459986161
That i have guessed at : 42.2075105, 20.7387093

Dataset extracted :
4673762 7478858 426
4674199 7478638 417
4673829 7478925 421
4674244 7479442 418
4673944 7479081 416
4673947 7479758 431
4674008 7479224 418
4673619 7478259 424
4673805 7478967 422
4673422 7472701 446
4673898 7479217 445
4673721 7478955 442
4674214 7479341 418
4673663 7477959 404
4674143 7478690 397
4672991 7475621 373
4672305 7472891 334
4672899 7480582 476
4674445 7479139 426
4674051 7479167 408

Attempts at swapping the xy and other things failed.

Command used :
cs2cs -E +proj=utm +zone=34T +units=m +proj=tmerc -f "%.7f"

example :
7478858 4673762 51.1551257 34.9253801 0.0000000
Updates of KPAOnline Scraper

here is the shell command : 

for x in resultPAmun.asp\?IS\=*PZ; do  perl ~/Desktop/maps/openstreetmapkosova/convert2.pl $x;  done > newtest2.osm

Updates of KPAOnline Scraper

here is my new script :

scraper: 

#street

#grep 'td width="189" bgcolor="#ECF5FF" class="fonteLatest"' $1;



########

# street 

#td width="189" bgcolor="#ECF5FF" class="fonteLatest"

#grep -h -A1 -e 'strong\>Street' $1

# | cut -d\> -f2 

grep -h -A1 -e 'Street' $1 | tail -1 | cut -d\> -f2 | cut -d\< -f1 > $1.street

 

## building (house number)

grep -h -A1 -e 'Building' $1 | tail -1 | cut -d\> -f2 | cut -d\< -f1 > $1.building



## GPS Grid UTM

#grep -h -A2 -e 'GPS' $1 | tail -1 | cut -d\> -f2 | cut -d\< -f1 > $1.gps

grep -h -A2 -e 'GPS' $1  | tail -1 |  cut -d\< -f1 | sed -e's;/; ;g' | cs2cs -E +proj=utm +zone=34T +units=m +proj=tmerc -f "%.7f"  > $1.gps



#perl ~/Desktop/maps/openstreetmapkosova/convert2osm.pl convert.txt > new.osm

converter: 

use strict;

use warnings;



my $filename=shift @ARGV;



sub getval

{

    my $self =shift;

    my $name =shift;

    open IN,$filename . ".$name" or die "no file";

    my $value;

#    warn "filename $filename";

    while ()

    {

#    warn $_;

	chomp;

	

	$value = $_;

	$value =~ s/\s+$//g;

    }

    close IN;

    $self->{attrs}{$name}=$value;

}



my $obj = {};



getval $obj, "street";

getval $obj, "building";

getval $obj, "gps";



my ($umtx,$umty,$lat,$lon)=split " ",$obj->{attrs}{gps};



$obj->{lat} = $lat;

$obj->{lon} = $lon;

$obj->{umtx} = $umtx;

$obj->{umty} = $umty;

delete $obj->{attrs}{gps};

if ($filename =~ /=(.+)/)

{

    $obj->{attrs}{kpaid}=$1;

}



#print q[



#my $id = -1;

my $oldnode = `grep $lat changeset.txt | grep lon`;



#warn $oldnode;



#print qq[

#

#];

#print ""



die "no old node" unless $oldnode =~ /node id/;



#$oldnode =~ s/node id/node chant/;

$oldnode =~ s/node id\=/node action="modify" id=/g;

$oldnode =~ s/\/\>/\>/g;

print $oldnode;



foreach my $k (sort keys %{$obj->{attrs}})

{

    my $v= $obj->{attrs}{$k};

    print "\n";

}



print "\n";



#print q[];

Updates of KPAOnline Scraper

Here is my changeset.
osm.org/browse/changeset/2077242
I removed one house from sharri. Hmm. Will have to do that another time.
nato attack sites in kosovo

Kosovo: Applying Geographic Information Systems in an International Humanitarian Crisis
http://proceedings.esri.com/library/userconf/proc00/professional/papers/PAP937/p937.htm

Kosovo: Applying GIS in an International Humanitarian Crisis
by David G. Smith, Department of State
http://www.esri.com/news/arcuser/0701/kosovo.html

Mineseeker airship :
http://maic.jmu.edu/Journal/5.1/Notes/Mineseeker/mineseeker.htm

Mine report kosovo :
http://www.mineaction.org/docs/1006648832_.asp

(open mine maps anyone)

Explosive ordnance disposal and mine clearance in Kosovo
http://www.reliefweb.int/rw/rwb.nsf/db900sid/MIRA-6EG484?OpenDocument

http://www.wsws.org/articles/2002/jun2002/trep-j28.shtml

http://www.washingtonpost.com/wp-dyn/content/article/2006/12/09/AR2006120900353_pf.html

http://www.icrc.org/web/eng/siteeng0.nsf/html/explosive-remnants-of-war-brochure-311201

http://articles.latimes.com/1999/jun/15/local/me-46797

'According to a KFOR map minefields are mainly present on the westerns slopes of the mountain range (confirming the report by jimorothy) but this map only depicts minefields which are known to UNMIK/KFOR and also it is not up to date (2005).
http://www.summitpost.org/mountain/rock/173204/-272-eravica.html#chapter_7

KFOR in mine awareness program
PRAGUE -- The Czech KFOR troops educate Kosovo children on how to avoid UXO and landmine injuries.

http://www.b92.net/eng/news/society-article.php?yyyy=2008&mm=10&dd=27&nav_id=54535

Yugoslavia: Education Programs Warn Kosovo's Children About Landmines
http://www.rferl.org/content/Article/1091705.html

GIS in the Kosovo ethnic conflict solution. The project "Sentinel".
http://proceedings.esri.com/library/userconf/proc00/professional/papers/PAP929/p929.htm

http://www.reliefweb.int/rw/rwb.nsf/db900SID/MHII-65952Q?OpenDocument
Military Technical Agreement: Between the International Security Force ("KFOR") and the Governments of the Federal Republic of Yugoslavia and the Republic of Serbia

http://www.defenselink.mil/news/newsarticle.aspx?id=42738
KFOR's First Priority: Countermine Operations

http://www.bulgaria-italia.com/fry/docs/Military_Agreement.htm
2. Phased Withdrawal of FRY Forces (ground): The FRY agrees to a phased withdrawal of all FRY Forces from Kosovo to locations in Serbia outside Kosovo. FRY Forces will mark and clear minefields, booby traps and obstacles. As they withdraw, FRY Forces will clear all lines of communication by removing all mines, demolitions, booby traps, obstacles and charges. They will also mark all sides of all minefields. International security forces' ("KFOR") entry and deployment into Kosovo will be synchronized. The phased withdrawal of FRY Forces from Kosovo will be in accordance with the sequence outlined below:

a. Detailed records, positions and descriptions of all mines, unexploded ordnance, explosive devices, demolitions, obstacles, booby traps, wire entanglement, physical or military hazards to the safe movement of any personnel in Kosovo laid by FRY Forces.

http://www.rta.nato.int/Pubs/RDP.asp?RDP=RTO-EN-SET-116
RTO-EN-SET-116 Low-Cost Navigation Sensors and Integration Technology

nato attack sites in kosovo

Lyx,
yes. There is also a huge number of mines left. good point!
screenscraping

here are my commands :
grep -h -A2 GPS resultPAmun.asp\?IS\=* | sort -u | grep ^04 | cut -d\< -f1 | sed -e 's;\/; ;g' > points.txt

cs2cs -E +proj=utm +zone=34T +units=m +proj=tmerc -f "%.9f" < points.txt > convert.txt
perl ~/Desktop/maps/openstreetmapkosova/convert2osm.pl convert.txt > new.osm

There were two data errors in the input :
0477588 4673224 20.798029599 4.281299370 0.000000000
0488275 4668363 20.868633114 36.760955805 0.000000000

my convert script : 

use strict;

use warnings;



print q[



  ];



my $id=0;

while (<>)

{



	if (/^\s*(\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)\s*$/)

	{	

	  $id--;

	  print qq[];

	}

else

{

die "error $_";

}



}



print q[];

screenscraping

here is my first point
osm.org/browse/changeset/2073369