#!/usr/bin/python ############################################################################ # # SCRIPT: ecw2tango.py Version 2 # Release two - some checks to see what units maps uses, we # do only degrees, not map coordinates. And degrees like # 148.0000000 nor 147d14'15.71"E, 36d40'29.14"S (buts in code) # # AUTHOR: David Bannon, Vic, Australia # Except two functions identified below obtained from OSM # # PURPOSE: Create map tiles (as per OSM) from a ECW map. # Needs gdal with ecw enabled. # # # COPYRIGHT: (c) 2009 David Bannon # Except two functions identified below obtained from OSM # # This program is free software under the GNU General Public # License (>=v2). Read the file publicly available GNU # documentation for further information. # # This script 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. Do not use # this script in any situation that could result in even the slightest # personal discomfort let alone risk to life or property. # ############################################################################# # # This short script takes an ECW format map and creates a set of map # tiles that are the equilivent (in format, not content) of those # produced by the OSM rendering engine. # I wrote it to make map tiles for off line use by TangoGPS but I # guess it may be useful for other purposes. While TangoGPS uses # OSM maps (and thats great!) we sometimes need alternatives - # * OSM maps consist mainly of roads, if you are going # where the roads are a bit thin, or want to see other # than roads, you need more. # * There are some very useful ECW maps published for # special places or applications. # # --- Use --- # ecw2tango [options] mapname.ecw # try ecw2tango --help # # --- Owners--- # The ECW format is owned by http://en.wikipedia.org/wiki/Leica_Geosystems see # http://www.erdas.com/Products/ERDASProductInformation/tabid/84/CurrentID/1142/Default.aspx # Its widely used for electronic maps, often survey or topographic. # The need to add ECW to gdal means you must build both from source, its easy. # http://trac.osgeo.org/gdal/wiki/ECW # Gdal is a translator library for raster geospatial data formats, see http://www.gdal.org/ # The OSM is a very cool open mapping project, http://www.osm.org # TangoGPS is a delightfully easy and very capable live map and GPS application # written by Marcus Bauer and friends, see http://www.tangogps.org # # LIMITATIONS # 1. This works only with maps in which the bitmap edges are parallel to the longitude # and lattitude lines. If the bitmap is bound by (eg) gridlines, its too hard (for me). # # # --- To Do --- # At some stage I might add two extra capabilities - # 1. Remove maps associated with a particular map. The ECW format seems heaps # kinder on disk space than PNG. # 2. Deal with boundries between maps. At present I start with first complete # tile to the right of the left edge etc. This means there will be a line of # missing tiles between each map. I reckon could write out an incomplete # tile, say 456R.png being the right hand side of tile 456 and if I see # and existing 456L.png, I'd stich them together. One day ... # ############################################################################### import commands import math import os import sys import string from optparse import OptionParser def GetToken(St, Token, First): # returns either the first or second number as a string # after the token. Always a comma between but any non numeric # is possible at the end. StepStart = St.find(Token) + len(Token) StepStop = S[StepStart:].find(",") if (First): XString = S[StepStart:StepStart+StepStop] return XString StepStart = StepStart + StepStop + 1 EndMarker = 0 Numerics = ['-',' ','.','0','1','2','3','4','5','6','7','8','9'] while S[StepStart+EndMarker] in Numerics: EndMarker = EndMarker + 1 XString = S[StepStart:StepStart+EndMarker] return XString # No, I don't understand. I don't have any maps that give dimensions # like this, so why did I put this function in ? # def XXConvertDegreeXX (St, Select): # 147d14'15.71"E, 36d40'29.14"S # 144d 0'0.00"E, 36d 0'0.00"S First = 1 while St[First] in string.digits: First = First + 1 Value = float(St[:First]) First = First + 1 Second = First Numerics = [' ','.','0','1','2','3','4','5','6','7','8','9'] while St[Second] in Numerics: Second = Second + 1 Value = Value + (float(St[First:Second]) / 60.0) Second = Second + 1 First = Second while St[Second] in Numerics: Second = Second + 1 Value = Value + (float(St[First:Second]) / 3600.0) Second = Second + 1 if (St[Second] == "W" or St[Second] == "S"): Value = -Value if Select == "Second": return ConvertDegree(St[Second + 2:], "First") return Value def ConvertDegree(St, Select): # 148.5000000, -37.0000000) End = 0 Numerics = ['-', ' ','.','0','1','2','3','4','5','6','7','8','9'] while St[End] in Numerics: End = End + 1 Value = float(St[:End]) if Select == "Second": return ConvertDegree(St[End+1:], "First") return Value def XXGetCornersXX(St, Search): SkipFirstBlock = 41 Start = St.find(Search) + SkipFirstBlock First = ConvertDegree(St[Start:Start+30], "First") if St[Start] == ",": Start = Start + 1 return First, ConvertDegree(St[Start:Start+30], "Second") def GetCorners(St, Search): SkipFirstBlock = 13 # thats where the data starts after the label. Start = St.find(Search) + SkipFirstBlock First = ConvertDegree(St[Start:Start+30], "First") if St[Start] == ",": Start = Start + 1 return First, ConvertDegree(St[Start:Start+30], "Second") # This function taken from OSM website # http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames # Give it a tile X and Y (in tile numb) and get back degrees def numb2deg(xtile, ytile, Zoom): n = 2.0 ** Zoom lon_deg = xtile / n * 360.0 - 180.0 lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n))) lat_deg = lat_rad * 180.0 / math.pi return (lon_deg, lat_deg) # This function taken from OSM website # http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames def deg2num(lon_deg, lat_deg, Zoom): lat_rad = lat_deg * math.pi /180.0 n = 2.0 ** Zoom xtile = int((lon_deg + 180.0) / 360.0 * n) ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n) return (xtile, ytile) # Writes out a 256x256 tile from a slice from the indicated map. def WritePNGFile(X, Y, TileName, XSlice, YSlice, MapName, Test): Xst = str(int(round(X))) Yst = str(int(round(Y))) CmdStr = "gdal_translate -srcwin " + Xst + " " + Yst + " " + str(XSlice) + " " + str(YSlice) + " -outsize 256 256 " + MapName + " -of png " + TileName + ".png" if options.Verbose: print "[", CmdStr, "]" if (not Test): if options.Verbose: Result = commands.getoutput(CmdStr) else: Result = commands.getoutput(CmdStr + " -quiet") if options.Verbose: print Result # Calculates the slice size out of the main map for a given zoom. # to make a 256x256 OSM tile. # Gets passed degrees but answers in pixels def SliceSize(XOff, YOff, Zoom): One = numb2deg(XOff, YOff, Zoom) Two = numb2deg(XOff+1, YOff+1, Zoom) DegXperTile = Two[0] - One[0] DegYperTile = One[1] - Two[1] fNoTilesX = (Right - Left)/DegXperTile fNoTilesY = (Top - Bottom)/DegYperTile X = int(round((float(XSize) / fNoTilesX))) Y = int(round((float(YSize) / fNoTilesY))) return (X, Y) # Give it a tile numbers and get back the top left pixel position # Uses Regional Vars, Right, Left, Top, Bottom (degrees), XSize, YSize (pixels), # and XOffset, YOffset (tiles) # XTile and YTile are relative into the map, start at 0,0 being # top, left (incomplete) tile. # If we are going to correct for maps that are not 'square', it might be here.. def SliceStart(XTile, YTile, Zoom, XTiles, YTiles): X, Y = numb2deg(XTile + XOffset, YTile + YOffset, Zoom) # ans in degrees XPixelPerDegree = XSize / (Right - Left) YPixelPerDegree = YSize / (Top - Bottom) return int(round((X - Left) * XPixelPerDegree)), int(round((Top - Y) * YPixelPerDegree)) # ----------------------------------------------------- # -------------- Start of main function --------------- # Right, lets get command line options first. # ----------------------------------------------------- parser = OptionParser() parser.add_option("-i", "--info", action="store_true", dest="Info", help="Just report on the map, don't write tiles") parser.add_option("-v", "--verbose", action="store_true", dest="Verbose", help="Show us all the noise as it happens.") parser.add_option("-t", "--test", action="store_true", dest="Test", help="Do a test run, tell us but dont write") parser.add_option("-d", "--dir_with_map_tile", action="store", type="string", dest="Dir", help="Where to write map tiles to") # parser.add_option("-r", "--remove_tiles", action=action="store_true", dest="Remove", # help="Remove tiles that relate to the given map.") parser.add_option("-z", "--zoom", action="store", type="int", dest="Zoom", help="Force a zoom level, 3-16") (options, args) = parser.parse_args() if len(args) == 0: print "No map file name provided, expected an ecw file" sys.exit(0) S=commands.getoutput("gdalinfo " + args[0]) UpperLeft = GetCorners(S, "Upper Left") UpperRight = GetCorners(S, "Upper Right") LowerLeft = GetCorners(S, "Lower Left") LowerRight = GetCorners(S, "Lower Right") # Some maps give those dimensions in degrees (ie 148.0000000), # others in map coordinates (ie 521239.115). Right now, we work # only with degrees. if abs(UpperLeft[0]) > 360: print "Sorry, cannot handle maps in map coordinates, only degrees" sys.exit(0) # OK, this is where we may fall down. # Identify if map is not 'square' wrt Lon, Lat # If the map is not bound by longitude and lattitude lines then I'm not # smart enough to do the necessary conversion ! if (int(round((UpperLeft[0] - LowerLeft[0]) * 1000.0)) <> 0) or \ (int(round((UpperRight[0] - LowerRight[0]) * 1000.0)) <> 0) or \ (int(round((UpperLeft[1] - UpperRight[1]) * 1000.0)) <> 0) or \ (int(round((LowerLeft[1] - LowerRight[1]) * 1000.0)) <> 0): print "Sorry, this does not look like a map we can convert." print "At this stage, all we can convert is maps that have" print "its boundries parallel to longitude and lattitude." sys.exit(0) XSize = int(GetToken(S, "Size is", True)) # answer in pixels. YSize = int(GetToken(S, "Size is", False)) Left = UpperLeft[0] # answer in degrees. Right = UpperRight[0] Top = UpperLeft[1] Bottom = LowerLeft[1] print "Map is ", XSize, "x", YSize print "Map Left", Left, " Right", Right, " Top", Top, " Bottom", Bottom print "X degrees per 256 pixels ", ((Left - Right)/XSize)*256 print "Y degress per 256 pixels ", ((Top - Bottom)/YSize)*256 Zoom = 0 if options.Zoom == None: Diff = 1000000 # ie just some big number ! # to guess zoom, we try and minimize difference between 256 and XSlice for Z in range(2, 16): Xs, Ys = SliceSize(0, 0, Z) if abs(256 - Xs) < Diff: Diff = abs(256 - Xs) Zoom = Z else: Zoom = int(options.Zoom) XOffset, YOffset = deg2num(Left, Top, Zoom) # answer in global tile numbers One = numb2deg(XOffset, YOffset, Zoom) Two = numb2deg(XOffset+1, YOffset+1, Zoom) XSlice, YSlice = SliceSize(XOffset, YOffset, Zoom) # answer is in pixels print "Choosen Zoom is ", Zoom, "Top left of first tile is", One print "Degrees per 256 pixel tile - lon", (One[0] - Two[0]), "lat", (One[1] - Two[1]) # print "XSlice is", XSlice, "YSlice is", YSlice, "number tiles", fNoTilesX, fNoTilesY print "XSlice is", XSlice, "YSlice is", YSlice if options.Dir == None: ZDir = str(Zoom) + "/" else: ZDir = options.Dir + "/" + str(Zoom) + "/" # OK, now we need the tile numbers we can generate from this map. Remember, we discard # any incomplete tiles for now. XTiles = 0 XTiles2 = 0 YTiles = 0 YTiles2 = 0 XPixels = 0 YPixels = 0 while XPixels < XSize: XTiles += 1 XPixels, YPixels = SliceStart(XTiles, 0, Zoom, 20, 20) XTiles -= 1 while YPixels < YSize: YTiles += 1 XPixels, YPixels = SliceStart(0, YTiles, Zoom, 20, 20) YTiles -= 1 YPixels = 0 while YPixels < YSize: YTiles2 += 1 XPixels, YPixels = SliceStart(XTiles, YTiles2, Zoom, XTiles, YTiles) YTiles -= 1 if YTiles2 < YTiles: YTiles = YTiles2 XPixels = 0 while XPixels < XSize: XTiles2 += 1 XPixels, YPixels = SliceStart(XTiles2, YTiles, Zoom, XTiles, YTiles) XTiles -= 1 if XTiles2 < XTiles: XTiles = XTiles2 print "Tile range will be ", XTiles, YTiles if options.Info: sys.exit(0) # OK, we now know all we need to run. # Create, if necessary, the Zoom directory try: # in case directory already exists. if (not options.Test): os.mkdir(ZDir) except OSError: if options.Verbose: print "directory exists", ZDir for Xt in range(1, XTiles): # X is column number, subdir to write into WDir = ZDir + str(XOffset + Xt) try: # in case directory already exists. if (not options.Test): os.mkdir(WDir) except OSError: if options.Verbose: print "Directory exists ", WDir print "Target Directory or Column = ", WDir, Xt, "of", XTiles for Yt in range(1, YTiles): # Y is file name if options.Verbose: print "We'd do ", WDir + "/" + str(Yt+YOffset) + ".png at", SliceStart(Xt, Yt, Zoom) X, Y = SliceStart(Xt, Yt, Zoom, XTiles, YTiles) WritePNGFile(X, Y, WDir + "/" + str(Yt + YOffset), XSlice, YSlice, args[0], options.Test) # Thats - X, Y pixels into bitmap, Full-file-name, size, size, input map, test-only print "done !" sys.exit(0) # Determine the location of the left edge of the second tile, we waste the incomplete first. # Similar for second down from the top. StartXPixels = ((XSlice) * (Two[0] - Left)) / (Two[0] - One[0]) StartYPixels = ((YSlice) * (Two[1] - Top)) / (Two[1] - One[1]) XOffset = XOffset + 1 YOffset = YOffset + 1 # OK, here we go. # XOffset and YOffset represent the OSM tile numbers (X is col, Y is file name) # for the choosen zoom. # XSlice and YSlice are the size of the bit we copy out on map. X = StartXPixels Y = StartYPixels # Create, if necessary, the Zoom directory try: # in case directory already exists. if (not options.Test): os.mkdir(ZDir) except OSError: if options.Verbose: print "directory exists", ZDir while X < float(XSize - XSlice): Y = StartYPixels YOff = YOffset WDir = ZDir + str(XOffset) try: # in case directory already exists. if (not options.Test): os.mkdir(WDir) except OSError: if options.Verbose: print "Directory exists ", WDir # if (not options.Test): os.chdir(str(XOffset)) print "Target Directory or Column = ", ZDir + str(XOffset) while Y < float(YSize - YSlice): WritePNGFile(X, Y, WDir + "/" + str(YOff), XSlice, YSlice, args[0], options.Test) YOff = YOff + 1 Y = Y + float(YSlice) # if (not options.Test): os.chdir("..") XOffset = XOffset + 1 X = X + float(XSlice) print "all done"