Merge branch 'develop'

Mostly fixes #1 although not elegant.

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 <input.tif> <corrected.tif>
+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] <kml>`
+
+Option | Result
+---------|--------
+`-h, --help` | Show this help message and exit
+`-o FILE, --outfile=FILE` | Write output to FILE
+`-v, --verbose` | Verbose output
+
+
+

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("""    <GroundOverlay>
                 <name>%s</name>
@@ -130,17 +171,18 @@ def create_kml(source, filename, directory, tile_size=1024, border=0, name=None,
     """ % bounds)
                 bob.write("""        </LatLonBox>
             </GroundOverlay>
-    """);
+    """)
 
     bob.write("""  </Folder>
     </kml>
     """)
 
     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)