extras/statistics/create-image.pl

#!/usr/bin/perl
#
# OSM data statistics Perl script
#
# Part of the Routino routing software.
#
# This file Copyright 2008-2016 Andrew M. Bishop
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

use Graphics::Magick;

# Range of tiles and zoom - PARAMETERS THAT CAN BE EDITED

@xrange=(120..129);             # At zoom level 8
@yrange=(73..87);               # At zoom level 8

$zbase=6;                       # Zoom level of images in basemap.

$z=13;                          # Granularity of data points.

# Base image dimensions

$basescale=2**(8-$zbase);

$width =(1+$#xrange)*256/$basescale;
$height=(1+$#yrange)*256/$basescale;

# Chosen zoom level

$scale=2**(8-$z);
$tilesize=256*$scale/$basescale;

# Get the command line arguments

$prefix=$ARGV[0];
$name  =$ARGV[1];
$number=$ARGV[2]+0;

# Graph annotations

$annotation{highway}  ="'highway=$name'";
$annotation{property} ="$name property";
$annotation{speed}    ="'maxspeed=$name'";
$annotation{transport}="$name allowed";

# Read in the data

%density=();

while(<STDIN>)
  {
   ($x,$y,@distances)=split(/ +/);

   $distance=$distances[$number];

   if($distance > 0)
     {
      $area=&xy_area($z,$x,$y);

      $density{$x,$y}=($distance/$area);
     }
  }

# Find the maximum value

$max=0;

foreach $xy (keys %density)
  {
   $density=$density{$xy};

   $max=$density if($density>$max);
  }

$max=500.0*int(($max+499)/500);
$max=500.0 if($max<500);

# Create the overlay image

$colour0=&colour(0);

$image=Graphics::Magick->new(size => "${width}x${height}");
$image->ReadImage("xc:$colour0");

foreach $xy (keys %density)
  {
   ($x,$y)=split($;,$xy);

   $colour=&colour($density{$x,$y}/$max);

   $x1=(($x*$scale)-$xrange[0])*256/$basescale;
   $y1=(($y*$scale)-$yrange[0])*256/$basescale;

   if($tilesize==1)
     {
      $image->Draw(primitive => 'point', points => "$x1,$y1",
                   fill => $colour,
                   antialias => 'false');
     }
   else
     {
      $x2=$x1+$tilesize-1;
      $y2=$y1+$tilesize-1;

      $image->Draw(primitive => 'rectangle', points => "$x1,$y1 $x2,$y2",
                   strokewidth => 0, stroke => $colour, fill => $colour,
                   antialias => 'false');
     }
  }

# Create the scale indicator

$indicator=Graphics::Magick->new(size => "${width}x40");
$indicator->ReadImage('xc:white');

foreach $v (0..($width-2*5))
  {
   $x=$v+5;
   $y1=12; $y2=23;
   $density=$v/($width-2*5);

   $indicator->Draw(primitive => 'line', points => "$x,$y1,$x,$y2",
                    stroke => &colour($density),
                    antialias => 'false');
  }

$indicator->Annotate(text => "0", font => 'Helvetica', pointsize => '12', style => Normal,
                     fill => 'black',
                     x => 5, y => 11, align => Left);

foreach $frac (0.25,0.5,0.75)
  {
   $indicator->Annotate(text => $max*$frac, font => 'Helvetica', pointsize => '12', style => Normal,
                        fill => 'black',
                        x => 5+($width-2*5)*$frac, y => 11, align => Center);
  }

$indicator->Annotate(text => $max, font => 'Helvetica', pointsize => '12', style => Normal,
                     fill => 'black',
                     x => $width-5, y => 11, align => Right);

$indicator->Annotate(text => "Highway density (metres/square km) for $annotation{$prefix} per zoom $z tile", font => 'Helvetica', pointsize => '14', style => Normal,
                     fill => 'black',
                     x => $width/2, y => 36, align => Center);

# Create the combined images

$base=Graphics::Magick->new;
$base->ReadImage("basemap.png");

$base->Composite(image => $image, compose => Dissolve, x => 0, y => 0, opacity => 50);

$final=Graphics::Magick->new(size => ($base->Get('width')+2)."x".($base->Get('height')+$indicator->Get('height')+3));
$final->ReadImage('xc:black');

$final->Composite(image => $base     , compose => Over, x => 1, y => 1                     );
$final->Composite(image => $indicator, compose => Over, x => 1, y => $base->Get('height')+2);

undef $base;
undef $image;
undef $indicator;

# Write out the images

print "Writing '$prefix-$name.png'\n";

$final->Write("$prefix-$name.png");

$final->Resize(width => $width/4, height => $final->Get('height')/4);
$final->Write("$prefix-$name-small.png");

undef $final;

#
# Subroutines
#

sub xy_ll_rad
{
 my($zoom,$x,$y)=@_;

 $PI=3.14159265358979323846;

 my($longitude)=$PI*(($x * (2**(1-$zoom)))-1);
 my($latitude)=($PI/2)*((4/$PI)*atan2(exp($PI*(1-($y * (2**(1-$zoom))))),1) - 1);

 return ($longitude,$latitude);
}


sub xy_area
{
 my($zoom,$x,$y)=@_;

 $RADIUS=6378.137;

 my($width,$height);

 if(defined $width{$y})
   {
    $width=$width{$y};
   }
 else
   {
    my($lon1,$lat1)=&xy_ll_rad($z,$x  ,$y);
    my($lon2,$lat2)=&xy_ll_rad($z,$x+1,$y);

    $width=$RADIUS*($lon2-$lon1)*cos($lat1);

    $width{$y}=$width;
   }

 if(defined $height{$y})
   {
    $height=$height{$y};
   }
 else
   {
    my($lon1,$lat1)=&xy_ll_rad($z,$x,$y  );
    my($lon2,$lat2)=&xy_ll_rad($z,$x,$y+1);

    $height=$RADIUS*($lat1-$lat2);

    $height{$y}=$height;
   }

 return $width*$height;
}

sub colour
{
 my($density)=@_;

 $density=sqrt($density);

 $density=0 if($density<0);
 $density=1 if($density>1);

 my($r,$g,$b);

 if($density<0.5)
   {
    $r=0;
    $g=int(255*(2*$density));
    $b=int(255*(1-2*$density));
   }
 else
   {
    $density-=0.5;

    $r=int(255*(2*$density));
    $g=int(255*(1-2*$density));
    $b=0;
   }

 sprintf("#%02x%02x%02x",$r,$g,$b);
}
1	#!/usr/bin/perl
2	#
3	# OSM data statistics Perl script
4	#
5	# Part of the Routino routing software.
6	#
7	# This file Copyright 2008-2016 Andrew M. Bishop
8	#
9	# This program is free software: you can redistribute it and/or modify
10	# it under the terms of the GNU Affero General Public License as published by
11	# the Free Software Foundation, either version 3 of the License, or
12	# (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 Affero General Public License for more details.
18	#
19	# You should have received a copy of the GNU Affero General Public License
20	# along with this program. If not, see <http://www.gnu.org/licenses/>.
21	#
22
23	use Graphics::Magick;
24
25	# Range of tiles and zoom - PARAMETERS THAT CAN BE EDITED
26
27	@xrange=(120..129); # At zoom level 8
28	@yrange=(73..87); # At zoom level 8
29
30	$zbase=6; # Zoom level of images in basemap.
31
32	$z=13; # Granularity of data points.
33
34	# Base image dimensions
35
36	$basescale=2**(8-$zbase);
37
38	$width =(1+$#xrange)*256/$basescale;
39	$height=(1+$#yrange)*256/$basescale;
40
41	# Chosen zoom level
42
43	$scale=2**(8-$z);
44	$tilesize=256*$scale/$basescale;
45
46	# Get the command line arguments
47
48	$prefix=$ARGV[0];
49	$name =$ARGV[1];
50	$number=$ARGV[2]+0;
51
52	# Graph annotations
53
54	$annotation{highway} ="'highway=$name'";
55	$annotation{property} ="$name property";
56	$annotation{speed} ="'maxspeed=$name'";
57	$annotation{transport}="$name allowed";
58
59	# Read in the data
60
61	%density=();
62
63	while(<STDIN>)
64	{
65	($x,$y,@distances)=split(/ +/);
66
67	$distance=$distances[$number];
68
69	if($distance > 0)
70	{
71	$area=&xy_area($z,$x,$y);
72
73	$density{$x,$y}=($distance/$area);
74	}
75	}
76
77	# Find the maximum value
78
79	$max=0;
80
81	foreach $xy (keys %density)
82	{
83	$density=$density{$xy};
84
85	$max=$density if($density>$max);
86	}
87
88	$max=500.0*int(($max+499)/500);
89	$max=500.0 if($max<500);
90
91	# Create the overlay image
92
93	$colour0=&colour(0);
94
95	$image=Graphics::Magick->new(size => "${width}x${height}");
96	$image->ReadImage("xc:$colour0");
97
98	foreach $xy (keys %density)
99	{
100	($x,$y)=split($;,$xy);
101
102	$colour=&colour($density{$x,$y}/$max);
103
104	$x1=(($x$scale)-$xrange[0])256/$basescale;
105	$y1=(($y$scale)-$yrange[0])256/$basescale;
106
107	if($tilesize==1)
108	{
109	$image->Draw(primitive => 'point', points => "$x1,$y1",
110	fill => $colour,
111	antialias => 'false');
112	}
113	else
114	{
115	$x2=$x1+$tilesize-1;
116	$y2=$y1+$tilesize-1;
117
118	$image->Draw(primitive => 'rectangle', points => "$x1,$y1 $x2,$y2",
119	strokewidth => 0, stroke => $colour, fill => $colour,
120	antialias => 'false');
121	}
122	}
123
124	# Create the scale indicator
125
126	$indicator=Graphics::Magick->new(size => "${width}x40");
127	$indicator->ReadImage('xc:white');
128
129	foreach $v (0..($width-2*5))
130	{
131	$x=$v+5;
132	$y1=12; $y2=23;
133	$density=$v/($width-2*5);
134
135	$indicator->Draw(primitive => 'line', points => "$x,$y1,$x,$y2",
136	stroke => &colour($density),
137	antialias => 'false');
138	}
139
140	$indicator->Annotate(text => "0", font => 'Helvetica', pointsize => '12', style => Normal,
141	fill => 'black',
142	x => 5, y => 11, align => Left);
143
144	foreach $frac (0.25,0.5,0.75)
145	{
146	$indicator->Annotate(text => $max*$frac, font => 'Helvetica', pointsize => '12', style => Normal,
147	fill => 'black',
148	x => 5+($width-25)$frac, y => 11, align => Center);
149	}
150
151	$indicator->Annotate(text => $max, font => 'Helvetica', pointsize => '12', style => Normal,
152	fill => 'black',
153	x => $width-5, y => 11, align => Right);
154
155	$indicator->Annotate(text => "Highway density (metres/square km) for $annotation{$prefix} per zoom $z tile", font => 'Helvetica', pointsize => '14', style => Normal,
156	fill => 'black',
157	x => $width/2, y => 36, align => Center);
158
159	# Create the combined images
160
161	$base=Graphics::Magick->new;
162	$base->ReadImage("basemap.png");
163
164	$base->Composite(image => $image, compose => Dissolve, x => 0, y => 0, opacity => 50);
165
166	$final=Graphics::Magick->new(size => ($base->Get('width')+2)."x".($base->Get('height')+$indicator->Get('height')+3));
167	$final->ReadImage('xc:black');
168
169	$final->Composite(image => $base , compose => Over, x => 1, y => 1 );
170	$final->Composite(image => $indicator, compose => Over, x => 1, y => $base->Get('height')+2);
171
172	undef $base;
173	undef $image;
174	undef $indicator;
175
176	# Write out the images
177
178	print "Writing '$prefix-$name.png'\n";
179
180	$final->Write("$prefix-$name.png");
181
182	$final->Resize(width => $width/4, height => $final->Get('height')/4);
183	$final->Write("$prefix-$name-small.png");
184
185	undef $final;
186
187	#
188	# Subroutines
189	#
190
191	sub xy_ll_rad
192	{
193	my($zoom,$x,$y)=@_;
194
195	$PI=3.14159265358979323846;
196
197	my($longitude)=$PI(($x (2**(1-$zoom)))-1);
198	my($latitude)=($PI/2)((4/$PI)atan2(exp($PI(1-($y (2**(1-$zoom))))),1) - 1);
199
200	return ($longitude,$latitude);
201	}
202
203
204	sub xy_area
205	{
206	my($zoom,$x,$y)=@_;
207
208	$RADIUS=6378.137;
209
210	my($width,$height);
211
212	if(defined $width{$y})
213	{
214	$width=$width{$y};
215	}
216	else
217	{
218	my($lon1,$lat1)=&xy_ll_rad($z,$x ,$y);
219	my($lon2,$lat2)=&xy_ll_rad($z,$x+1,$y);
220
221	$width=$RADIUS($lon2-$lon1)cos($lat1);
222
223	$width{$y}=$width;
224	}
225
226	if(defined $height{$y})
227	{
228	$height=$height{$y};
229	}
230	else
231	{
232	my($lon1,$lat1)=&xy_ll_rad($z,$x,$y );
233	my($lon2,$lat2)=&xy_ll_rad($z,$x,$y+1);
234
235	$height=$RADIUS*($lat1-$lat2);
236
237	$height{$y}=$height;
238	}
239
240	return $width*$height;
241	}
242
243	sub colour
244	{
245	my($density)=@_;
246
247	$density=sqrt($density);
248
249	$density=0 if($density<0);
250	$density=1 if($density>1);
251
252	my($r,$g,$b);
253
254	if($density<0.5)
255	{
256	$r=0;
257	$g=int(255(2$density));
258	$b=int(255(1-2$density));
259	}
260	else
261	{
262	$density-=0.5;
263
264	$r=int(255(2$density));
265	$g=int(255(1-2$density));
266	$b=0;
267	}
268
269	sprintf("#%02x%02x%02x",$r,$g,$b);
270	}