Index: src/map_download.c =================================================================== --- src/map_download.c (revision 2476) +++ src/map_download.c (revision 2477) @@ -61,8 +61,8 @@ # endif -#include -#include +#include "gpsdrive.h" +#include "map_handler.h" #include "gpsdrive_config.h" #include "gui.h" #include "map_download.h" @@ -133,24 +133,20 @@ MAPSOURCE_LANDSAT, "1 : 50 million", 0, 50000000, MAPSOURCE_OSM_TAH, "OpenStreetMap Tiles@Home", -1, -1, - MAPSOURCE_OSM_TAH, "1 : 147 456 000", 1, 256*576000, - MAPSOURCE_OSM_TAH, "1 : 73 728 000", 2, 128*576000, - MAPSOURCE_OSM_TAH, "1 : 36 864 000", 3, 64*576000, - MAPSOURCE_OSM_TAH, "1 : 18 432 000", 4, 32*576000, - MAPSOURCE_OSM_TAH, "1 : 9 216 000", 5, 16*576000, - MAPSOURCE_OSM_TAH, "1 : 4 608 000", 6, 8*576000, - MAPSOURCE_OSM_TAH, "1 : 2 304 000", 7, 4*576000, - MAPSOURCE_OSM_TAH, "1 : 1 152 000", 8, 2*576000, - - MAPSOURCE_OSM_TAH, "1 : 576 000", 9, 576000, - MAPSOURCE_OSM_TAH, "1 : 288 000", 10, 288000, - MAPSOURCE_OSM_TAH, "1 : 144 000", 11, 144000, - MAPSOURCE_OSM_TAH, "1 : 72 000", 12, 72000, - MAPSOURCE_OSM_TAH, "1 : 36 000", 13, 36000, - MAPSOURCE_OSM_TAH, "1 : 18 000", 14, 18000, - MAPSOURCE_OSM_TAH, "1 : 9 000", 15, 9000, - MAPSOURCE_OSM_TAH, "1 : 4 500", 16, 4500, - MAPSOURCE_OSM_TAH, "1 : 2 250", 17, 2250, + /* scale varies with latitude, so this is just a rough guide + which will only be valid for mid-lats */ + /* Octave code: for lat=0:5:75; disp( [lat (a * 2*pi *pixelfact * cos(lat * pi/180)) / (256*2^9)]); end */ + MAPSOURCE_OSM_TAH, "1 : 2 500", 17, 2250, + MAPSOURCE_OSM_TAH, "1 : 5 000", 16, 4500, + MAPSOURCE_OSM_TAH, "1 : 10 000", 15, 9000, + MAPSOURCE_OSM_TAH, "1 : 20 000", 14, 18000, + MAPSOURCE_OSM_TAH, "1 : 40 000", 13, 36000, + MAPSOURCE_OSM_TAH, "1 : 75 000", 12, 72000, + MAPSOURCE_OSM_TAH, "1 : 150 000", 11, 144000, + MAPSOURCE_OSM_TAH, "1 : 300 000", 10, 288000, + MAPSOURCE_OSM_TAH, "1 : 600 000", 9, 576000, + /* the distortion gets too bad at scales wider than 1:500k + It would be more accurate to stop earlier, but we compromise */ MAPSOURCE_N_ITEMS, "", 0, 0 }; @@ -165,7 +161,23 @@ gtk_entry_set_text (GTK_ENTRY (lon_entry), lon); } +/* ***************************************************************************** + * calculate the local map scale based on Web Tile zoom level and latitude + */ +double +calc_webtile_scale (double lat, int zoom) +{ + double scale; + double a = 6378137.0; /* major radius of WGS84 ellipsoid (meters) */ + scale = (a * 2*M_PI * cos(DEG2RAD(lat)) * PIXELFACT) / (256 * pow(2,zoom)); + + if (mydebug > 3) + g_print ("Tile scale: %.2f\n", scale); + + return scale; +} + /* ***************************************************************************** * callback to set paramaters for map to download */ @@ -195,9 +207,13 @@ 3, &mapdl_scale, -1); /* TODO: determine map_ or top_ proj at this point so drawdownloadrectangle() * knows how big to draw the green preview box */ + if (local_config.mapsource_type == MAPSOURCE_OSM_TAH) + mapdl_scale = (int)calc_webtile_scale(mapdl_lat, mapdl_zoom); + if (mydebug > 3) g_print ("new map scale/zoom level: %d / %d\n", mapdl_scale, mapdl_zoom); + local_config.mapsource_scale = gtk_combo_box_get_active (GTK_COMBO_BOX (scale_combobox)); current.needtosave = TRUE; } @@ -439,6 +455,7 @@ gchar file_path[512]; gchar path[40]; gchar img_fmt[4]; + gchar scale_str[40]; if (mapdl_active) return TRUE; @@ -450,11 +467,14 @@ g_strlcpy (path, "landsat", sizeof (path)); g_strlcpy (img_fmt, "jpg", sizeof (img_fmt)); mapdl_geturl_landsat (); + g_snprintf (scale_str, sizeof (scale_str), "%d", mapdl_scale); break; case MAPSOURCE_OSM_TAH: g_strlcpy (path, "openstreetmap_tah", sizeof (path)); g_strlcpy (img_fmt, "png", sizeof (img_fmt)); mapdl_geturl_osm_tah (); + mapdl_scale = (int)calc_webtile_scale(mapdl_lat, mapdl_zoom); + g_snprintf (scale_str, sizeof (scale_str), "%d", mapdl_zoom); break; default: return TRUE; @@ -464,10 +484,11 @@ g_print (" download url:\n%s\n", mapdl_url); /* set file path and create directory if necessary */ - - - g_snprintf (file_path, sizeof (file_path), "%s%s/%d/%.0f/%.0f/", - local_config.dir_maps, path, mapdl_scale, mapdl_lat, mapdl_lon); + + g_snprintf (file_path, sizeof (file_path), "%s%s/%s/%.0f/%.0f/", + local_config.dir_maps, path, scale_str, mapdl_lat, mapdl_lon); + + if(!g_file_test (file_path, G_FILE_TEST_IS_DIR)) { if (g_mkdir_with_parents (file_path, 0700)) @@ -477,14 +498,15 @@ } /* complete filename */ - g_snprintf (mapdl_file_w_path, sizeof (mapdl_file_w_path), "%s%s_%d_%5.3f_%5.3f.%s", - file_path, mapdl_proj, mapdl_scale, mapdl_lat, mapdl_lon, img_fmt); + g_snprintf (mapdl_file_w_path, sizeof (mapdl_file_w_path), "%s%s_%s_%5.3f_%5.3f.%s", + file_path, mapdl_proj, scale_str, mapdl_lat, mapdl_lon, img_fmt); if (mydebug > 10) g_print (" filename: %s\n", mapdl_file_w_path); gtk_progress_bar_set_text (GTK_PROGRESS_BAR (mapdl_progress), _("Loading Map...")); + mapdl_abort = FALSE; mapdl_download (); Index: src/map_download.h =================================================================== --- src/map_download.h (revision 2476) +++ src/map_download.h (revision 2477) @@ -1,10 +1,9 @@ /*********************************************************************** Copyright (c) 2008 Guenther Meyer - Website: www.gpsdrive.de/ -Disclaimer: Please do not use for navigation. +Disclaimer: Please do not use as a primary means of navigation. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -36,5 +35,6 @@ void mapdl_set_coords (gchar *lat, gchar *lon); gint mapdl_init (void); gint mapdl_cleanup (void); +double calc_webtile_scale (double, int); #endif /* GPSDRIVE_MAP_DOWNLOAD_H */ Index: scripts/gpsfetchmap.pl =================================================================== --- scripts/gpsfetchmap.pl (revision 2479) +++ scripts/gpsfetchmap.pl (revision 2480) @@ -47,6 +47,7 @@ modified (Jul 2007) by Maciek Kaliszewski modified (Oct 2007) by Andreas Putzo modified (Jan 2008) by Gernot Hillier (added Openstreetmap support) +modified (Aug 2009) by Hamish Bowman (fix OSM+WMS map scales) Version svn-$Version "; @@ -151,24 +152,20 @@ 10000000 =>10000000, 50000000 =>50000000, }, +# Web Map Tiles scale varies with latitude by the formula +# (a * 2*pi * PIXELFACT * cos(lat * pi/180)) / (256 * 2^zoom) +# where a is major radius of the Earth according to the WGS84 ellipsoid definition. +# zoom levels lower than 9 (~1:500k) badly distort the Mercator so are not offered. openstreetmap_tah => { - 256*576000 => 1, - 128*576000 => 2, - 64*576000 => 3, - 32*576000 => 4, - 16*576000 => 5, - 8*576000 => 6, - 4*576000 => 7, - 2*576000 => 8, - 576000 => 9, - 288000 => 10, - 144000 => 11, - 72000 => 12, - 36000 => 13, - 18000 => 14, - 9000 => 15, - 4500 => 16, - 2250 => 17 + 600000 => 9, + 300000 => 10, + 150000 => 11, + 75000 => 12, + 40000 => 13, + 20000 => 14, + 10000 => 15, + 5000 => 16, + 2500 => 17 } }; @@ -240,8 +237,6 @@ $polite=1; } - - if ( $mapserver eq 'geoscience' ) { $scale ||= join(",",keys %{$Scale2Zoom->{'geoscience'}}); @@ -269,22 +264,24 @@ pod2usage(-verbose=>2) if $man; sub append_koords($$$$); # {} +sub calc_webtile_scale ($$); # {} sub check_coverage($); # {} sub check_koord_file($); # {} sub debug($); # {} +sub file_count($); # {} sub geoscience_url($$$); # {} -sub landsat_url($$$); # {} -sub resize($$); #{} -sub file_count($); # {} sub get_coords_for_route; # {} sub get_coords_for_track($); # {} +sub get_gpsd_position(); # {} sub get_waypoint($); # {} sub is_map_file($); # {} +sub landsat_url($$$); # {} +sub openstreetmap_tah_url($$$); # {} sub read_gpstool_map_file(); # {} sub read_koord_file($); # {} +sub resize($$); # {} sub update_gpsdrive_map_koord_file(); # {} sub wget_map($$$); # {} -sub get_gpsd_position(); # {} STDERR->autoflush(1); STDOUT->autoflush(1); @@ -407,7 +404,9 @@ my ($existing,$wanted) = file_count($desired_locations); print "You are about to download $wanted (".($existing+$wanted).") file(s).\n"; - +if ($debug) { + print "Politeness delay set to $polite seconds,\n"; +} if ( $mapserver eq 'geoscience' ){ print "+-----------------------------------------------------------+\n"; print "| Geoscience Maps are Copyright, Commonwealth of Australia |\n"; @@ -419,6 +418,7 @@ print "+----------------------------------------------------------------+\n"; print "| Landsat maps are courtesy JPL/NASA's OnEarth WMS Global Mosaic |\n"; print "| Map prefix will be set automatically based on scale. |\n"; + print "+----------------------------------------------------------------+\n"; # By law, US Government data is without copyright. }elsif ( $mapserver eq 'openstreetmap_tah' ){ print "+-----------------------------------------------------------+\n"; @@ -427,6 +427,8 @@ print "| They are free for use under the terms of the |\n"; print "| Creative Commons \"Attribution-Share Alike 2.0 Generic\" |\n"; print "| license. See http://www.openstreetmap.org for details. |\n"; + print "+-----------------------------------------------------------+\n"; + } elsif ( ! $force) { print "You may violating the map servers terms of use!\n"; print "Please use their service responsible!\n"; @@ -654,14 +656,32 @@ ###################################################################### sub map_filename($$$){ my ($scale,$lati,$long) = @_; + + my $mscale = $scale; + + if ( $mapserver eq 'openstreetmap_tah' ) { + my $zoom = undef; + for my $s ( sort keys %{$Scale2Zoom->{openstreetmap_tah}} ) { + next unless $s == $scale; + $zoom = $Scale2Zoom->{openstreetmap_tah}->{$s}; + last; + } + unless ( $zoom ) { + print "Error calculating Zoomlevel for Scale: $scale\n"; + return (undef); + } + $mscale = $zoom; + } - my $filename = "$mapserver/$scale" + my $filename = "$mapserver/$mscale" ."/".int($lati) # ."/".sprintf("%3.1f",$lati) ."/".int($long) - ."/$FILEPREFIX$scale-$lati-$long.$fileext"; - printf("Filename(%.0f,%.5f,%.5f): $filename\n",$scale,$lati,$long) + ."/$FILEPREFIX$mscale-$lati-$long.$fileext"; + + printf("Filename(%.0f,%.5f,%.5f): $filename\n",$mscale,$lati,$long) if $debug; + return $filename; } @@ -678,6 +698,9 @@ # mapblast/1000/047/047.0232/9/map_1000-047.0232-0009.8140.gif 47.02320 9.81400 1000 #my $filename = "$mapserver/$scale/".int($lati)."/".sprintf("%3.1f",$lati). #"/".int($long)."/$FILEPREFIX$scale-$lati-$long.gif"; + if ($debug) { + print "---------------\n"; + } # redundant?? if ($mapserver eq 'landsat') @@ -1082,7 +1105,7 @@ $lon2 += 360; } - debug( "landsat_url(LAT=$lat,LON=$lon,SCALE=$scale,FACTOR=$factor)"); + debug( "landsat_url(LAT=$lat,LON=$lon,SCALE=$scale)"); debug( "Calculated Lat1 $lat1"); debug( "Calculated Lat2 $lat2"); @@ -1111,7 +1134,28 @@ return ($url,$scale); } + ############################################################################# +# calculate the local map scale based on Web Tile zoom level and latitude +# This same formula is used by OSM, Google Maps, Microsoft Maps, ... +sub calc_webtile_scale ($$){ + my $lat = shift; + my $zoom = shift; + + # major radius of WGS84 ellipsoid (meters) + my $a = $RADIUS_KM * 1000; + + my $zscale = (2*$a*$PI * cos($lat*$D2R) * $PIXELFACT) / (256 * 2**$zoom); + + if ($debug) { + print "Tile scale: $zscale\n"; + } + + return ($zscale); +} + + +############################################################################# sub openstreetmap_tah_url($$$){ my $lati = shift; my $long = shift; @@ -1122,7 +1166,7 @@ for my $s ( sort keys %{$Scale2Zoom->{openstreetmap_tah}} ) { next unless $s == $scale; $zoom = $Scale2Zoom->{openstreetmap_tah}->{$s}; - $mapscale = $s; + $mapscale = calc_webtile_scale($lati, $zoom); last; } @@ -1132,10 +1176,9 @@ } if ($debug) { - print "\n"; - print "Using openstreetmap_tah zoom ", $zoom, " for requested scale ", $scale, ":1 actual scale ", $mapscale, ":1\n"; - print "lat: $lati\n"; - print "lon: $long\n"; + print "Using openstreetmap_tah zoom ", $zoom, " for requested scale 1:", $scale, ". Actual scale 1:", $mapscale, "\n"; + print "cetner lat: $lati\n"; + print "cetner lon: $long\n"; } my $url = "http://tah.openstreetmap.org/MapOf/?lat=$lati&long=$long&z=$zoom&w=1280&h=1024&format=png"; @@ -1903,7 +1946,7 @@ gpsfetchmap -w -sc -a <#> -p -gpsfetchmap -la -lo -sc -a <#> -p +gpsfetchmap -la -lo -sc -a <#> -p gpsfetchmap -sla -endla -slo -endlo -sc -a <#> -p @@ -1915,7 +1958,7 @@ [-la ] [-lo ] [-sla ] [-endla ] [-slo ] [-endlo ] - [-sc ] [-a <#>] [-p] [-m ] + [-sc ] [-a <#>] [-p <#>] [-m ] [-u ] [-md ] [-W ] [-t ] [-r] [-C ] [-P ] [-F] [-d] [-v] [-h] [-M] [-n] [-U] [-c] @@ -2000,7 +2043,7 @@ of longitude. 'units' is read from the configuration file (-C) or as defined by (-u). -=item B<-p, --polite> +=item B<-p, --polite <#>> This causes the program to sleep one second between downloads to be polite to the mapserver. Takes an optional value of number of seconds to sleep. @@ -2117,6 +2160,13 @@ Prints the manual page and exits. +=back + + +=head1 OUTPUT + +=over 2 + =item B When downloading Maps the output reads as folows: @@ -2130,5 +2180,16 @@ =back + +=head1 EXAMPLE + +Download all 1:10000 OpenStreetMaps in a 5 km radius around the +Sydney Convention Centre. Between tile downloads it will let +the server rest for 3 seconds. + +gpsfetchmap --mapserver openstreetmap_tah --polite 3 \ + --lat "-33.8753" --lon 151.2001 --scale 10000 \ + --area 5 --unit kilometers + =cut Index: scripts/mapnik/gpsdrive_mapnik_gentiles-in.py =================================================================== --- scripts/mapnik/gpsdrive_mapnik_gentiles-in.py (revision 2473) +++ scripts/mapnik/gpsdrive_mapnik_gentiles-in.py (revision 2479) @@ -7,14 +7,20 @@ Options: -h, --help show this help - -b, --bbox boundingbox (lon,lat,lon,lat) - Be carefull! Quote negative values! - -s, --scale scale single/range - 1 = 147456000 - 2 = 73728000 - ... - 15 = 9000 - 16 = 4500 - 17 = 2250 + -b, --bbox boundingbox (minlon,minlat,maxlon,maxlat) + - Be carefull! "Quote" negative values! + -s, --scale scale single/range (zoom level or min-max) + (below "9" Mercator becomes distorted; + actual scale will vary with latitude) + 9 - 1:600,000 + 10 - 1:300,000 + 11 - 1:150,000 + 12 - 1:75,000 + 13 - 1:40,000 + 14 - 1:20,000 + 15 - 1:10,000 + 16 - 1:5,000 + 17 - 1:2,500 --test testrun = generates Munich example Examples: @@ -22,7 +28,7 @@ Munich: gpsdrive_mapnik_gentiles.py -b 11.4,48.07,11.7,48.2 -s 10-16 - World: + World: (just to demonstrate lat/lon order; scale not recommended) gpsdrive_mapnik_gentiles.py -b "-180.0,-90.0,180.0,90.0" -s 1-6 """ @@ -32,12 +38,19 @@ import getopt import string -zoom2scale = [728 * 576000,256 * 576000,128 * 576000,64 * 576000,32 * 576000,16 * 576000,8 * 576000,4 * 576000,2 * 576000,576000,288000,144000,72000,36000,18000,9000,4500,2250,1125] - DEG_TO_RAD = pi/180 RAD_TO_DEG = 180/pi +def calc_scale (lat, zoom): + # GpsDrive's hardcoded pixels per meter ratio + PixelFact = 2817.947378 + # wgs84 major Earth axis + a = 6378137.0 + dynscale = ( a * 2*pi * cos(lat * DEG_TO_RAD) * PixelFact ) / ( 256*pow(2,zoom) ) + print dynscale + return dynscale + def minmax (a,b,c): a = max(a,b) a = min(a,c) @@ -90,17 +103,28 @@ os.mkdir(tile_dir) gprj = GoogleProjection(maxZoom+1) + #m = Map(2 * 256,2 * 256) m = Map(1280,1024) load_map(m,mapfile) - prj = Projection("+proj=merc +datum=WGS84") - + + #prj = Projection("+proj=merc +datum=WGS84") + # What follows is from /usr/share/proj/esri.extra, EPSG:900913 + # "Chris' funny epsgish code for the google mercator" + prj = Projection("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs") + ll0 = (bbox[0],bbox[3]) ll1 = (bbox[2],bbox[1]) - + for z in range(minZoom,maxZoom + 1): + if z == 9: + print "CAUTION: Mercator projection begins to be noticeably distorted at this zoom level." + elif z < 9: + print "WARNING: Mercator projection is very distorted at this zoom level." + px0 = gprj.fromLLtoPixel(ll0,z) px1 = gprj.fromLLtoPixel(ll1,z) + for x in range(int(px0[0]/640.0),int(px1[0]/640.0)+1): for y in range(int(px0[1]/512.0),int(px1[1]/512.0)+1): p0 = gprj.fromPixelToLL((x * 640.0, (y+1) * 512.0),z) @@ -145,7 +169,7 @@ fh_mapkoord.write(tile_path + " ") fh_mapkoord.write(str((p0[1] + p1[1]) / 2) + " ") fh_mapkoord.write(str((p0[0] + p1[0]) / 2) + " ") - fh_mapkoord.write(str(zoom2scale[z])) + fh_mapkoord.write(str( calc_scale( (p0[1] + p1[1])/2, z) )) fh_mapkoord.write(" " + str(p0[1]) + " " + str(p0[0])) fh_mapkoord.write(" " + str(p1[1]) + " " + str(p1[0])) fh_mapkoord.write("\n") @@ -221,6 +245,7 @@ # check for correct values if str(eval(bboxs[0])) != bboxs[0] or str(eval(bboxs[1])) != bboxs[1] or str(eval(bboxs[2])) != bboxs[2] or str(eval(bboxs[3])) != bboxs[3]: + # rounding problems... what exactly is this supposed to be checking ??? sys.exit("Boundingbox invalid!") if minZoom < 1 or minZoom > 17 or maxZoom < 1 and maxZoom > 17 or minZoom > maxZoom or int(minZoom) <> minZoom or int(maxZoom) <> maxZoom: