source: subversion/applications/utils/osm-extract/polygons/extract-polygon_0.4.pl

Last change on this file was 5344, checked in by frederik, 11 years ago

0.4 versions not used anymore

  • Property svn:executable set to *
File size: 6.8 KB
Line 
1#!/usr/bin/perl
2
3# This script extracts all nodes, segments, and ways lying inside a polygon
4#
5# The script first compiles a polygon from the data, then simplifies it,
6# and then processes the planet file. In processing the file, first a
7# bounding box is used to exclude all nodes that lie outside, and the
8# more expensive polygon check is only made for those inside.
9#
10# Author Frederik Ramm <frederik@remote.org> / public domain
11# Author Keith Sharp <kms@passback.co.uk> / public domain
12
13# Progress meter, exclusions, et al all added by Martijn van Oosterhout
14# <kleptog@svana.org> for the AND import (August 2007) / public domain
15
16# Get this from CPAN!
17use Math::Polygon;
18use Getopt::Long;
19use File::Basename;
20use strict;
21
22my $borderpolys = [];
23my $currentpoints;
24my $minlon = 999;
25my $minlat = 999;
26my $maxlat = -999;
27my $maxlon = -999;
28
29my $used_nodes = {};
30my $used_segments = {};
31
32my $help;
33my $polyfile;
34my $infile;
35my $outfile;
36my $remainsfile;
37my $compress;
38my $verbose = 0;
39
40my $prog = basename($0);
41
42sub usage ()
43{
44        print STDERR << "EOF";
45
46This Perl script will process a planet.osm file and extract the nodes,
47segments, and ways, falling within a polygon.
48
49usage: $prog [-h] [-c number] [-i file] [-o file] [-r file] -p file
50
51 -h       : print ths help message and exit.
52 -c num   : reduce the number of points in the polygon by num %.  This 
53            only takes effect if the polgon has more 100 points.
54 -p file  : file containing the polygon definition.
55 -i file  : OSM planet file to process.
56 -o file  : OSM planet file to output.
57 -r file  : (optional) outputs a simple list of all node ids that are within
58            your selected area but are also used by segments not selected.
59            Some forms of post-processing need this information.
60
61If you do not supply an input or output file input will be read from STDIN
62and output written to STDOUT.
63
64A polygon file should be a plain text file with each line containing a
65longitude followed by a latitude, for example a polygon defining Great
66Britain and Ireland (and some of the smaller islands) would be:
67
681
69        -0.6450E+01     0.4980E+02
70        -0.2000E+01     0.4890E+02
71        -0.1850E+01     0.4925E+02
72        -0.2080E+01     0.4973E+02
73         0.1350E+01     0.5090E+02
74         0.2250E+01     0.5258E+02
75        -0.0500E+01     0.6130E+02
76        -0.8920E+01     0.5785E+02
77        -0.1140E+02     0.5130E+02
78        -0.6450E+01     0.4980E+02
79END
80
81Multiple polygons can be described by adding additional sections to the
82polygon file.
83
84EOF
85        exit;
86}
87
88GetOptions ('h|help' => \$help,
89                        'c|compress=i' => \$compress,
90                        'v|verbose' => \$verbose,
91                        'i|input=s' => \$infile,
92                        'o|output=s' => \$outfile,
93                        'p|polygon=s' => \$polyfile,
94                        'r|remains=s' => \$remainsfile) || usage ();
95
96if ($help) {
97        usage ();
98}
99
100if (! $polyfile) {
101        usage ();
102}
103
104open (PF, "$polyfile") || die "Could not open file: $polyfile: $!";
105
106my $invert;
107# initialize border polygon.
108while(<PF>)
109{
110    if (/^(!?)\d/)
111    {
112        $invert = ($1 eq "!") ? 1 : 0;
113        $currentpoints = [];
114    }
115    elsif (/^END/)
116    {
117        my $pol = Math::Polygon->new(points => $currentpoints);
118        if (($compress > 0 && $compress < 100) && $pol->nrPoints > 99) {
119                my $simp = $pol->simplify($pol->nrPoints*(100-$compress)/100);
120                push(@{$borderpolys}, [$simp,$invert,[$simp->bbox]]);
121        } else {
122                push(@{$borderpolys}, [$pol,$invert,[$pol->bbox]]);
123                }
124        undef $currentpoints;
125        if( $verbose )
126        {
127            printf STDERR "Added polygon: %d points (%.2f,%.2f)-(%.2f,%.2f) %s\n",
128                $borderpolys->[-1][0]->nrPoints,
129                @{$borderpolys->[-1][2]},
130                ($borderpolys->[-1][1] ? "exclude" : "include");
131        }
132    }
133    elsif (defined($currentpoints))
134    {
135        /^\s+([0-9.E+-]+)\s+([0-9.E+-]+)/ or die;
136        push(@{$currentpoints}, [$2, $1]);
137        $minlat = $2 if ($2 < $minlat);
138        $maxlat = $2 if ($2 > $maxlat);
139        $minlon = $1 if ($1 < $minlon);
140        $maxlon = $1 if ($1 > $maxlon);
141    }
142}
143
144close (PF);
145
146if ($outfile) {
147        open (OF, ">$outfile") || die "Could not open file: $outfile: $!";
148} else {
149        *OF = *STDOUT;
150}
151
152if ($infile) {
153        if ($infile =~ /\.bz2/ ) {
154                open (IF, "bzcat $infile |") || die "Could not bzcat file: $infile: $!";
155        } else {
156                open (IF, "<$infile") || die "Could not open file: $infile: $!";
157        }
158} else {
159        *IF = *STDIN;
160}
161
162if ($remainsfile) {
163        open (RF, ">$remainsfile") || die "Could not open file: $remainsfile: $!";
164}
165
166print OF << "EOF";
167<?xml version="1.0" encoding="UTF-8"?>
168<osm version="0.3" generator="extract-polygon.pl">
169EOF
170
171my $copy;
172my $waybuf;
173my $count = 0;
174
175while(<IF>) 
176{
177    $count++;
178    if( $verbose and ($count % 10000) == 0 )
179    {
180        # When reading from a compressed file we assume a compression ration of 14.6
181        my $perc = tell(IF)*100/((-s IF) || ((-s $infile)*14.6));
182        if( $perc > 100 ) { $perc = 100 }
183        printf STDERR "\r%.2f%% ", $perc;
184    }
185    last if /^\s*<\/osm>/;
186
187    # Note: we allow a minus in the ID to allow incomplete files also   
188    if (/^\s*<node.*id=["'](-?\d+)['"].*lat=["']([Ee0-9.-]+)["'] lon=["']([Ee0-9.-]+)["']/)
189    {
190        $copy = 0;
191
192        next if (($2 < $minlat) || ($2 > $maxlat));
193        next if (($3 < $minlon) || ($3 > $maxlon));
194
195        my $ll = [$2, $3];
196
197        foreach my $p (@{$borderpolys})
198        {
199            my($poly,$invert,$bbox) = @$p;
200            next if ($ll->[0] < $bbox->[0]) or ($ll->[0] > $bbox->[2]);
201            next if ($ll->[1] < $bbox->[1]) or ($ll->[1] > $bbox->[3]);
202
203            if ($poly->contains($ll))
204            {
205                # If this polygon is for exclusion, we immediately bail and go for the next point
206                if( $invert )
207                {
208                    $copy = 0;
209                    last;
210                }
211                $copy = 1;
212            }
213        }
214        if( $copy )
215        {
216            $used_nodes->{$1} = 1;
217            print OF;
218        }
219    }
220    elsif (/^\s*<segment id=['"](-?\d+)["'] from=['"](-?\d+)["'] to=["'](-?\d+)["']/)
221    {
222        $copy = 0;
223        if ($used_nodes->{$2} && $used_nodes->{$3})
224        {
225            $used_segments->{$1}=1;
226            print OF;
227            $copy = 1;
228        }
229        elsif ($remainsfile)  # Only test if we want the output
230        {
231            if ($used_nodes->{$2})
232            {
233                print RF "$2\n";
234            }
235            if ($used_nodes->{$3})
236            {
237                print RF "$3\n";
238            }
239        }
240    }
241    elsif (/^\s*<way /)
242    {
243        $copy = 0;
244        $waybuf = $_;
245        undef $used_nodes;
246    }
247    elsif ($copy)
248    {
249        print OF;
250    }
251    elsif (/^\s*<seg id=['"](.*)["']/)
252    {
253        if ($used_segments->{$1})
254        {
255            print OF $waybuf;
256            print OF;
257            $copy = 1;
258        }
259        else
260        {
261            $waybuf .= $_;
262        }
263    }
264}
265
266print OF << "EOF";
267</osm>
268EOF
269
270close (OF);
271close (IF);
272close (RF) if ($remainsfile);
273print STDERR "\n" if $verbose;
274
275exit;
Note: See TracBrowser for help on using the repository browser.