diff --git a/README.md b/README.md new file mode 100644 index 0000000..c3f815d --- /dev/null +++ b/README.md @@ -0,0 +1,75 @@ +## GDAL to CustomMap Tools +**NOTE**: *The majority of this code comes from [`gdal2custommap`](https://github.com/tf198/gdal2custommap) by [Tris Forster](https://github.com/tf198) written for Python 2. This fork is made to work with Python 3 and also fixes some bugs with large raster maps.* + +Some python scripts to generate Garmin CustomMap compatible KMZ from +geo-referenced images. + +Garmin CustomMap has some restrictions on the KMZ it can display + +1. No tile larger than 1024x1024 (or equivalent total pixels) +2. No more than 100 tiles on the device (e.g. 4 KMZ files with 25 tiles each) + +`gdal2tiles` (the GDAL default tiler) uses 256x256 tiles so requires 4 times +as many tiles as necessary. `gdal2kml` tries to tile the map as efficiently +as possible so you dont end up with thin tiles at the right or bottom. + +`gdal2kml` **DOES NOT** do any warping so the source image must be correctly georeferenced and in `EPSG:4326` +projection or the conversion will fail. +For other projections you need to get the GDAL utility programs and run it through `gdal_warp` first. If you are not sure about the projection run it through anyway as it will not adversly affect the output. +This may introduce a black border which you can then cut out using the --crop option to gdal2kml. +``` +gdalwarp -t_srs EPSG:4326 -r cubic +gdal2kml.py corrected.tif output.kml --crop 20 +``` + +### Requirements +You need the GDAL library and Python bindings installed. On Ubuntu +its as simple as +``` +sudo apt-get install python-gdal +``` + +For windows there are prebuilt binaries at http://www.gisinternals.com +~~The current stable is 1.9.0 and I have tested with the MSVC2008 version and the +following packages but it can be a bit tricky to get all the paths right~~ +``` +(http://python.org/download) +python-2.7.2.msi +(http://www.gisinternals.com/sdk/PackageList.aspx?file=release-1500-gdal-1-9-mapserver-6-0.zip) +gdal-19-1500-core.msi +gdal-19-1500-ecw.msi +GDAL-1.9.0.win32-py2.7.msi +``` + +### Basic Usage +```bash +python gdal2kml.py input.tif my_map.kml +python kml2kmz.py my_map.kml +``` + +### gdal2kml.py + +Usage: `gdal2kml.py [options] src_file dst_file` + +Option | Result +---------|-------- +`-h, --help` | show this help message and exit +`-d DIRECTORY, --dir=DIRECTORY` | Where to create jpeg tiles +`-c BORDER, --crop=BORDER` | Crop border +`-n NAME, --name=NAME` | KML folder name for output +`-o ORDER, --draw-order=ORDER` | KML draw order +`-t TILE_SIZE, --tile-size=TILE_SIZE` | Max tile size [1024] +`-q QUALITY, --quality=QUALITY` | JPEG output quality 0-100 [75] +`-v, --verbose` | Verbose output + +### kml2kmz.py +Usage: `kml2kmz.py [options] ` + +Option | Result +---------|-------- +`-h, --help` | Show this help message and exit +`-o FILE, --outfile=FILE` | Write output to FILE +`-v, --verbose` | Verbose output + + + diff --git a/gdal2kml.py b/gdal2kml.py index c108ebe..24cff36 100755 --- a/gdal2kml.py +++ b/gdal2kml.py @@ -1,91 +1,124 @@ -#!/usr/bin/env python - -import os, math, re, sys, subprocess, logging +import os +import math +import logging from optparse import OptionParser -from osgeo import gdal +from osgeo import gdal +from osgeo import osr + def tiles(canvas, target=1024): - """ + """ Brute force algorithm to determine the most efficient tiling method for a given canvas If anyone can figure out a prettier one please let me know - is actually harder then you'd think! """ best_case = (canvas[0] * canvas[1]) / float(target**2) - + # handle the trivial cases first - if canvas[0] <= target: return [ 1, math.ceil(best_case) ] - if canvas[1] <= target: return [ math.ceil(best_case), 1 ] - - r = [ float(x) / target for x in canvas ] - - # brute force the 4 methods - - a_up = [ math.ceil(x) for x in r ] - b_up = [ math.ceil(best_case / x) for x in a_up ] - - a_down = [ math.floor(x) for x in r ] - b_down = [ math.ceil(best_case / x) for x in a_down ] - + if canvas[0] <= target: + return [1, math.ceil(best_case)] + + if canvas[1] <= target: + return [math.ceil(best_case), 1] + + r = [float(x) / target for x in canvas] + + # brute force the 4 methods + + a_up = [math.ceil(x) for x in r] + b_up = [math.ceil(best_case / x) for x in a_up] + + a_down = [math.floor(x) for x in r] + b_down = [math.ceil(best_case / x) for x in a_down] + results = [] for i in range(2): results.append((a_up[i], b_up[i], a_up[i] * b_up[i])) results.append((a_down[i], b_down[i], a_down[i] * b_down[i])) - + results.sort(key=lambda x: x[2]) - return [ int(x) for x in results[0][0:2] ] + return [int(x) for x in results[0][0:2]] -def create_tile(source, filename, offset, size, quality=75): + +def transform(x, y, geotransform): + xt = geotransform[0] + x*geotransform[1] + y*geotransform[2] + yt = geotransform[3] + x*geotransform[4] + y*geotransform[5] + + return xt, yt + + +def create_tile(img, filename, offset, size, quality=75): """ Create a jpeg of the given area and return the bounds. """ mem_drv = gdal.GetDriverByName('MEM') - mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) - bands = range(1, source.RasterCount+1) + mem_ds = mem_drv.Create('', size[0], size[1], img.RasterCount) + bands = range(1, img.RasterCount + 1) - data = source.ReadRaster(offset[0], offset[1], size[0], size[1], size[0], size[1], band_list=bands) + # TODO: consider this instead https://rasterio.readthedocs.io/en/latest/ + data = img.ReadRaster(offset[0], offset[1], size[0], size[1], size[0], size[1], band_list=bands) mem_ds.WriteRaster(0, 0, size[0], size[1], data, band_list=bands) + # Error comes because we go out of bounds of image? jpeg_drv = gdal.GetDriverByName('JPEG') jpeg_ds = jpeg_drv.CreateCopy(filename, mem_ds, strict=0, options=["QUALITY={0}".format(quality)]) - t = source.GetGeoTransform() - if t[2]!=0 or t[4]!=0: raise Exception("Source projection not compatible") - def transform((x, y)): - return ( t[0] + x*t[1] + y*t[2], t[3] + x*t[4] + y*t[5] ) - - nw = transform(offset) - se = transform([ offset[0] + size[0], offset[1] + size[1] ]) - + geotransform = img.GetGeoTransform() + + if geotransform[2] != 0 or geotransform[4] != 0: + raise Exception('Source projection incompatible, transform contains rotation') + else: + nw = transform(offset[0], offset[1], geotransform) + se = transform(offset[0] + size[0], offset[1] + size[1], geotransform) + result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } - + jpeg_ds = None mem_ds = None - + return result -def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): + +def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=None, quality=75): """ Create a kml file and associated images for the given georeferenced image """ + if exclude is None: + exclude = [] + img = gdal.Open(source) - img_size = [img.RasterXSize, img.RasterYSize] + projection = img.GetProjection() + logging.info(projection) - logging.debug('Image size: %s' % img_size) + # https://gdal.org/user/raster_data_model.html#raster-data-model + srs = osr.SpatialReference(wkt=projection) + authority = (srs.GetAttrValue('AUTHORITY', 0), srs.GetAttrValue('AUTHORITY', 1)) + logging.debug(authority) + + if authority != ('EPSG', '4326'): + # https://gdal.org/tutorials/osr_api_tut.html#coordinate-transformation + # ct = osr.CoordinateTransformation() + errmsg = 'Input file is not in standard CRS. Should be EPSG 4326 but is %s %s' % (authority[0], authority[1]) + logging.error(errmsg) + raise NotImplementedError(errmsg) - cropped_size = [ x - border * 2 for x in img_size ] + img_size = [img.RasterXSize, img.RasterYSize] + logging.debug('Image size: %s' % img_size) + cropped_size = [x - border * 2 for x in img_size] base, ext = os.path.splitext(os.path.basename(source)) - if not name: name = base + if not name: + name = base path = os.path.relpath(directory, os.path.dirname(filename)) tile_layout = tiles(cropped_size, tile_size) - tile_sizes = [ int(math.ceil(x)) for x in [ cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1] ] ] + tile_sizes = [int(math.floor(x)) for x in [cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1]]] logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') @@ -104,12 +137,20 @@ def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, logging.debug("Excluding tile %s" % tile) else: src_corner = (border + t_x * tile_sizes[0], border + t_y * tile_sizes[1]) - src_size = [ tile_sizes[0], tile_sizes[1] ] - if src_corner[0] + tile_sizes[0] > img_size[0] - options.border: src_size[0] = int(tile_sizes[0]) - if src_corner[1] + tile_sizes[1] > img_size[1] - options.border: src_size[1] = int(tile_sizes[1]) + src_size = [tile_sizes[0], tile_sizes[1]] + if src_corner[0] + tile_sizes[0] > img_size[0] - border: + src_size[0] = int(tile_sizes[0]) + + if src_corner[1] + tile_sizes[1] > img_size[1] - border: + src_size[1] = int(tile_sizes[1]) - outfile = "%s_%d_%d.jpg" % (base, t_x, t_y) - bounds = create_tile(img, "%s/%s" % (directory, outfile), src_corner, src_size, quality) + outfile = '%s_%d_%d.jpg' % (base, t_x, t_y) + outpath = '%s/%s' % (directory, outfile) + + if src_corner[0] + src_size[0] > img_size[0]: + logging.error('Pixel range outside image data!') + logging.error('Image width %i, trying to get at x=%i' % (img_size[0], src_corner[0] + src_size[0])) + bounds = create_tile(img, outpath, src_corner, src_size, quality) bob.write(""" %s @@ -130,7 +171,7 @@ def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, """ % bounds) bob.write(""" - """); + """) bob.write(""" @@ -138,9 +179,10 @@ def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, bob.close() img = None - -if __name__=='__main__': - usage = "usage: %prog [options] src_file dst_file" + + +if __name__ == '__main__': + usage = 'usage: %prog [options] src_file dst_file' parser = OptionParser(usage) parser.add_option('-d', '--dir', dest='directory', help='Where to create jpeg tiles') parser.add_option('-c', '--crop', default=0, dest='border', type='int', help='Crop border') @@ -152,28 +194,41 @@ def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, options, args = parser.parse_args() - if len(args) != 2: parser.error('Missing file paths') - source, dest = args + if len(args) != 2: + parser.error('Missing file paths') + + source_file, destination_file = args - if options.verbose: logging.basicConfig(level=logging.DEBUG) + if options.verbose: + logging.basicConfig(level=logging.DEBUG) # validate a few options - if not os.path.exists(source): parser.error('unable to file src_file') - #if options.scale<10 or options.scale>150: parser.error('scale must be between 10% and 150%') + if not os.path.exists(source_file): + parser.error('unable to file src_file') + # if options.scale<10 or options.scale>150: parser.error('scale must be between 10% and 150%') # set the default folder for jpegs - if not options.directory: options.directory = "%s.files" % os.path.splitext(dest)[0] - if not os.path.exists(options.directory): os.mkdir(options.directory) + if not options.directory: + options.directory = "%s.files" % os.path.splitext(destination_file)[0] + + if not os.path.exists(options.directory): + os.mkdir(options.directory) + logging.info('Writing jpegs to %s' % options.directory) # load the exclude file - exclude_file = source + ".exclude" + exclude_file = source_file + ".exclude" exclude = [] if os.path.exists(exclude_file): logging.debug("Using exclude file %s" % exclude_file) for line in open(exclude_file): exclude.append(line.rstrip()) logging.debug(exclude) - - create_kml(source, dest, options.directory, - tile_size=options.tile_size, border=options.border, name=options.name, order=options.order, exclude=exclude, quality=options.quality) + + create_kml(source_file, destination_file, options.directory, + tile_size=options.tile_size, + border=options.border, + name=options.name, + order=options.order, + exclude=exclude, + quality=options.quality) diff --git a/readme.rst b/readme.rst deleted file mode 100644 index 674dccd..0000000 --- a/readme.rst +++ /dev/null @@ -1,79 +0,0 @@ -GDAL to CustomMap Tools -======================= -Some python scripts to generate Garmin CustomMap compatible KMZ from -geo-referenced images. - -Garmin CustomMap has some restrictions on the KMZ it can display - -1. No tile larger than 1024x1024 (or equivalent total pixels) -2. No more than 100 tiles on the device (e.g. 4 KMZ files with 25 tiles each) - -``gdal2tiles`` (the GDAL default tiler) uses 256x256 tiles so requires 4 times -as many tiles as necessary. ``gdal2kml`` tries to tile the map as efficiently -as possible so you dont end up with thin tiles at the right or bottom. - -``gdal2kml`` **DOES NOT** do any warping so the source image must be correctly georeferenced and in ``EPSG:4326`` -projection or the conversion will fail. -For other projections you need to get the GDAL utility programs and run it through ``gdal_warp`` first. If you -are not sure about the projection run it through anyway as it will not adversly affect the output. -This may introduce a black border which you can then cut out using the --crop option to gdal2kml. -:: - - gdalwarp -t_srs EPSG:4326 -r cubic - gdal2kml.py corrected.tif output.kml --crop 20 - -Requirements ------------- -You need the GDAL library and Python bindings installed. On Ubuntu -its as simple as:: - - sudo apt-get install python-gdal - -For windows there are prebuilt binaries at http://www.gisinternals.com -The current stable is 1.9.0 and I have tested with the MSVC2008 version and the -following packages but it can be a bit tricky to get all the paths right:: - - (http://python.org/download) - python-2.7.2.msi - (http://www.gisinternals.com/sdk/PackageList.aspx?file=release-1500-gdal-1-9-mapserver-6-0.zip) - gdal-19-1500-core.msi - gdal-19-1500-ecw.msi - GDAL-1.9.0.win32-py2.7.msi - -Basic Usage ------------ -:: - - python gdal2kml.py input.tif my_map.kml - python kml2kmz.py my_map.kml - -gdal2kml.py ------------ -Usage: gdal2kml.py [options] src_file dst_file - -Options: - -h, --help show this help message and exit - -d DIRECTORY, --dir=DIRECTORY - Where to create jpeg tiles - -c BORDER, --crop=BORDER - Crop border - -n NAME, --name=NAME KML folder name for output - -o ORDER, --draw-order=ORDER - KML draw order - -t TILE_SIZE, --tile-size=TILE_SIZE - Max tile size [1024] - -q QUALITY, --quality=QUALITY JPEG output quality 0-100 [75] - -v, --verbose Verbose output - -kml2kmz.py ----------- -Usage: kml2kmz.py [options] - -Options: - -h, --help show this help message and exit - -o FILE, --outfile=FILE - Write output to FILE - -v, --verbose Verbose output - - -