source: subversion/applications/utils/ooc/reproject.pl @ 29291

Last change on this file since 29291 was 20278, checked in by richard, 10 years ago

harder better faster stronger

File size: 5.0 KB
Line 
1#!/usr/bin/perl -w
2
3        use Geo::Proj4;
4        use Math::Trig;
5        use Imager;
6        use Imager::CountColor;
7        use POSIX qw(floor ceil);
8
9        $google=Geo::Proj4->new(proj=>"merc", a=>6378137, b=>6378137,
10                                                    lat_ts=>0.0, lon_0=>0.0, x_0=>0.0, y_0=>0,
11                                                    k=>1.0, units=>"m", nadgrids=>'@null');
12        $blank=Imager::Color->new(0,255,0);
13        $scale1=15;     $scale2=16;     # Target zoom levels (min/max)
14        $root="tiles/";
15
16        # Do we need to reproject any images?
17        @files=();
18        for $fn (@ARGV) {
19                if ($fn=~/_p/) {
20                        # already reprojected, so no
21                        push @files,$fn;
22                } else {
23                        print "Reprojecting $fn\n";
24                        $fn2=$fn; $fn2=~s/_r/_p/; push @files,$fn2;
25                        system "gdalwarp --config GDAL_CACHEMAX 500 -wm 500 -s_srs EPSG:27700 -t_srs EPSG:900913 -srcnodata '0 255 0' -dstnodata '0 255 0' $fn $fn2";
26                        unlink $fn;
27                }
28        }
29
30        # Combine images
31        print "Combining\n";
32        system "gdalbuildvrt -srcnodata '0 255 0' combined.vrt @files";
33        system "gdal_translate combined.vrt combined.tiff";
34
35        # Get bounds
36        $info=`gdalinfo combined.tiff`;
37        if ($info=~/Size is (\d+), (\d+)/i) { $width=$1; $height=$2; }
38        if ($info=~/Upper Left\s*\(\s*([\d\-.]+),\s*([\d\-.]+)/i) { $xleft=$1; $ytop=$2; }
39        if ($info=~/Lower Right\s*\(\s*([\d\-.]+),\s*([\d\-.]+)/i) { $xright=$1; $ybottom=$2; }
40        if ($info=~/Pixel Size = \(-?([\d.]+)/i) { $pixelsize=$1; }
41       
42        ($umaxlat,$uminlong)=$google->inverse($xleft ,$ytop   );
43        ($uminlat,$umaxlong)=$google->inverse($xright,$ybottom);
44        print "Mercator: ($xleft, $ytop) - ($xright, $ybottom)\n";
45        print "Lat/long: ($umaxlat, $uminlong) - ($uminlat, $umaxlong)\n";
46
47        for ($scale=$scale1; $scale<=$scale2; $scale++) {
48                print "\nScale $scale\n";
49
50                # Get tiles
51                $gx1=long2tile($uminlong)+1; $gx2=long2tile($umaxlong)-1;
52                $gy1=lat2tile ($umaxlat )+1; $gy2=lat2tile ($uminlat )-1;
53                $gwidth =$gx2-$gx1; $gwidthpx =$gwidth *256;
54                $gheight=$gy2-$gy1; $gheightpx=$gheight*256;
55       
56                # Convert back to lat/long so we can crop
57                $gminlong=tile2long($gx1); $gmaxlong=tile2long($gx2);
58                $gmaxlat =tile2lat ($gy1); $gminlat =tile2lat ($gy2);
59       
60                print "\nGoogle crop range:\n";
61                print "$gminlong, $gmaxlong\n";
62                print "$gminlat, $gmaxlat\n";
63
64                print "\nGoogle tile range:\n";
65                print "X: $gx1 to $gx2 ($gwidth, $gwidthpx pixels)\n";
66                print "Y: $gy1 to $gy2 ($gheight, $gheightpx pixels)\n";
67
68                ($gprojleft ,$gprojtop   )=$google->forward($gmaxlat,$gminlong);
69                ($gprojright,$gprojbottom)=$google->forward($gminlat,$gmaxlong);
70                print "\nProjected Google values:\n";
71                print "Mercator: ($gprojleft, $gprojtop) - ($gprojright, $gprojbottom)\n";
72
73                $crop_left  =($gprojleft-$xleft)/$pixelsize;
74                $crop_right =($xright-$gprojright)/$pixelsize;
75                $crop_top   =($ytop-$gprojtop)/$pixelsize;
76                $crop_bottom=($gprojbottom-$ybottom)/$pixelsize;
77                print "crop: $crop_left, $crop_right, $crop_top, $crop_bottom\n";
78       
79                $ewidth=$width-$crop_left-$crop_right;
80                $eheight=$height-$crop_top-$crop_bottom;
81                print "expected width: $ewidth\n";
82                print "expected height: $eheight\n";
83       
84                # ------------------------
85                # Create tiles with Imager
86
87                # Crop to relevant size
88                # bbox of original is $uminlong->$umaxlong, $uminlat->$umaxlat
89                # what we want is     $gminlong->$gmaxlong, $gminlat->$gmaxlat
90
91                $pxsize=$width/($umaxlong-$uminlong);   # pixels per degree
92                print "\nResize:\n";
93                print "pxsize: $pxsize\n";
94
95                $crop_left  =($gminlong-$uminlong)*$pxsize;
96                $crop_right =($umaxlong-$gmaxlong)*$pxsize;
97                $crop_top   =(lat2y($umaxlat)-lat2y($gmaxlat))*$pxsize;
98                $crop_bottom=(lat2y($gminlat)-lat2y($uminlat))*$pxsize;
99
100                $owidth=$ewidth/$gwidthpx*256;
101                $oheight=$eheight/$gheightpx*256;
102                print "Source tile width $owidth, source tile height $oheight\n";
103
104                print "\nProcessing file:\n";
105                $pr=Imager->new();
106                $pr->open(file=>"combined.tiff") or die "Open error: ".$pr->errstr;
107
108                print "crop: $crop_left, $crop_right, $crop_top, $crop_bottom\n";
109                $pr=$pr->crop(left  =>$crop_left,
110                                          right =>$width-$crop_right,
111                                          top   =>$crop_top,
112                                          bottom=>$height-$crop_bottom) or die "Crop error: ".$pr->errstr;
113
114                print "write\n";
115                for ($x=$gx1; $x<$gx2; $x++) {
116                        for ($y=$gy1; $y<$gy2; $y++) {
117                                unless (-d "${root}$scale") { mkdir "${root}$scale"; }
118                                unless (-d "${root}$scale/$x") { mkdir "${root}$scale/$x"; }
119                                $tile=$pr->crop(left=>($x-$gx1)*$owidth , width=>$owidth,
120                                                                top =>($y-$gy1)*$oheight, height=>$oheight) or die "Tile crop error: ".$pr->errstr;
121                                $tile=$tile->scale(xpixels=>256, ypixels=>256, type=>'nonprop') or die "Tile scale error: ".$pr->errstr;
122                                $bcount=Imager::CountColor::count_color($tile,$blank);
123                                if ($bcount<5000) {
124                                        $tile->write(file=>"${root}$scale/$x/$y.jpg", jpegquality=>90) or die "Tile write error: ".$tile->errstr;
125                                }
126                        }
127                }
128        }
129        unlink "combined.tiff";
130        unlink "combined.vrt";
131       
132        print chr(7)."Finished\n";
133
134        # ==========
135        # Projection
136
137        sub long2tile { return (floor(($_[0]+180)/360*2**$scale)); }
138        sub lat2tile { return (floor((1-log(tan($_[0]*pi/180) + 1/cos($_[0]*pi/180))/pi)/2 *2**$scale)); }
139        sub tile2long { return ($_[0]/2**$scale*360-180); }
140        sub tile2lat { my $n=pi-2*pi*$_[0]/2**$scale; return (180/pi*atan(0.5*(exp($n)-exp(-$n)))); }
141        sub lat2y { return 180/pi*log(tan(pi/4+$_[0]*(pi/180)/2)); }
Note: See TracBrowser for help on using the repository browser.