source: subversion/applications/rendering/mapnik/all_tiles/tilecount.pl @ 34619

Last change on this file since 34619 was 1713, checked in by jonb, 13 years ago

Fix tilecount.pl calculation for zoom 17 & 18

  • Property svn:executable set to *
File size: 2.8 KB
Line 
1#!/usr/bin/perl
2use strict;
3#-----------------------------------------------------------------------------
4# says which tiles contain data from a specified OSM XML file
5# usage: tilecount.pl planet.osm > tile_list.txt
6#-----------------------------------------------------------------------------
7# Copyright 2006, Oliver White
8#
9# This program is free software; you can redistribute it and/or
10# modify it under the terms of the GNU General Public License
11# as published by the Free Software Foundation; either version 2
12# of the License, or (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program; if not, write to the Free Software
21# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
22#-----------------------------------------------------------------------------
23use Math::Trig;
24# Setup coordinate system
25my $LimitY = ProjectF(85.0511);
26my $LimitY2 = ProjectF(-85.0511);
27my $RangeY = $LimitY - $LimitY2;
28
29my %Tiles;
30
31# precalculate number of tiles in each zoom-level
32my @Num;
33foreach my $Z(0..18){
34  $Num[$Z] = 2 ** $Z;
35}
36   
37# Parse the planet.osm file
38open IN, "<",shift() || die("Need an OSM file");
39while(my $Line = <IN>){
40  if($Line =~ /<node.*lat=\"(.*?)\" lon=\"(.*?)\"/){
41    #print "$1, $2\n";
42    HandlePos($1,$2);
43  }
44}
45close IN;
46
47# Display the results
48foreach my $Key (keys %Tiles){
49  print "$Key\n";
50}
51
52sub HandlePos(){
53  my ($Lat,$Long) = @_;
54  # Find position in mercator projection
55  my $MercY = ProjectF($Lat);
56#  my $RelY = ($MercY - $LimitY2) / $RangeY;
57  my $RelY = ( $LimitY - $MercY ) / $RangeY;
58  my $RelX = ($Long + 180) / 360;
59 
60  # Convert to tile coordinates
61  foreach my $Z (0..18){
62    my $X = int($RelX * $Num[$Z]);
63    my $Y = int($RelY * $Num[$Z]);
64    my $ID = "$Z:$X:$Y";
65    $Tiles{$ID} = 1;
66    }
67}
68
69sub ProjectF($){
70  my $Lat = DegToRad(shift());
71  my $Y = log(tan($Lat) + sec($Lat));
72  return($Y);
73}
74sub Project(){
75  my ($Y, $Zoom) = @_;
76 
77  my $Unit = 1 / (2 ** $Zoom);
78  my $relY1 = $Y * $Unit;
79  my $relY2 = $relY1 + $Unit;
80 
81  $relY1 = $LimitY - $RangeY * $relY1;
82  $relY2 = $LimitY - $RangeY * $relY2;
83   
84  my $Lat1 = ProjectMercToLat($relY1);
85  my $Lat2 = ProjectMercToLat($relY2);
86  return(($Lat1, $Lat2)); 
87}
88sub ProjectMercToLat($){
89  my $MercY = shift();
90  return(RadToDeg(atan(sinh($MercY))));
91}
92sub ProjectL(){
93  my ($X, $Zoom) = @_;
94 
95  my $Unit = 360 / (2 ** $Zoom);
96  my $Long1 = -180 + $X * $Unit;
97  return(($Long1, $Long1 + $Unit)); 
98}
99sub DegToRad($){return pi * shift() / 180;}
100sub RadToDeg($){return 180 * shift() / pi;}
101
Note: See TracBrowser for help on using the repository browser.