Skip to main content
deleted 10569 characters in body
Source Link
Buff Fox
  • 305
  • 1
  • 14

A more manual method is the below python code that I adapted from https://github.com/tf198/gdal2custommap

https://github.com/Mp3Robbo/GeoTIFF_To_Garmin_KMZ

This uses PyQGIS to make a .kmz with tilesWhich I adapted from

import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'ReduRes' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 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[0] = tile_sizes[0] - 1 tile_sizes[1] = tile_sizes[1] - 1 print('Tile size is') print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name) for t_y in range(tile_layout[1]): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) if tile in exclude: 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] - 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) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile) 

https://github.com/tf198/gdal2custommap

Though OkMap does seem easier, and I'm goingThis uses PyQGIS to give itmake a go...kmz with tiles

A more manual method is the below code that I adapted from https://github.com/tf198/gdal2custommap

This uses PyQGIS to make a .kmz with tiles

import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'ReduRes' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 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[0] = tile_sizes[0] - 1 tile_sizes[1] = tile_sizes[1] - 1 print('Tile size is') print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name) for t_y in range(tile_layout[1]): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) if tile in exclude: 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] - 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) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile) 

Though OkMap does seem easier, and I'm going to give it a go...

A more manual method is the below python code

https://github.com/Mp3Robbo/GeoTIFF_To_Garmin_KMZ

Which I adapted from

https://github.com/tf198/gdal2custommap

This uses PyQGIS to make a .kmz with tiles

added 132 characters in body
Source Link
Buff Fox
  • 305
  • 1
  • 14
import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'BigBoyMap''ReduRes' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 path = os.path.relpath(directory, os.path.dirname(filename)) tile_layout = tiles(cropped_size, tile_size)  #swapsie = tile_layout[0] #tile_layout[0] = tile_layout[1] #tile_layout[1] = swapsie tile_sizes = [ int(math.ceil(x)) for x in [ cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1] ] ]  tile_sizes[0] = tile_sizes[0] - 1 tile_sizes[1] = tile_sizes[1] - 1   print('Tile size is')  print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name)   for t_y in range(tile_layout[1]-1): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) print(tile) if tile in exclude: 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] - 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)  print (src_corner) bounds = create_tile(img, "%s/%s" % (directory, outfile), src_corner, src_size, quality) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile)   
import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'BigBoyMap' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 path = os.path.relpath(directory, os.path.dirname(filename)) tile_layout = tiles(cropped_size, tile_size)  #swapsie = tile_layout[0] #tile_layout[0] = tile_layout[1] #tile_layout[1] = swapsie tile_sizes = [ int(math.ceil(x)) for x in [ cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1] ] ] print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name) for t_y in range(tile_layout[1]-1): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) print(tile) if tile in exclude: 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] - 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)  print (src_corner) bounds = create_tile(img, "%s/%s" % (directory, outfile), src_corner, src_size, quality) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile)   
import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'ReduRes' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 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[0] = tile_sizes[0] - 1 tile_sizes[1] = tile_sizes[1] - 1   print('Tile size is')  print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name)   for t_y in range(tile_layout[1]): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) if tile in exclude: 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] - 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) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile) 
deleted 24 characters in body
Source Link
Buff Fox
  • 305
  • 1
  • 14
import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'BigBoyMap' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':directory + finalName + '.tif'source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 path = os.path.relpath(directory, os.path.dirname(filename)) tile_layout = tiles(cropped_size, tile_size) #swapsie = tile_layout[0] #tile_layout[0] = tile_layout[1] #tile_layout[1] = swapsie tile_sizes = [ int(math.ceil(x)) for x in [ cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1] ] ] print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name) for t_y in range(tile_layout[1]-1): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) print(tile) if tile in exclude: 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] - 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) print (src_corner) bounds = create_tile(img, "%s/%s" % (directory, outfile), src_corner, src_size, quality) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile) 
import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'BigBoyMap' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':directory + finalName + '.tif'}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 path = os.path.relpath(directory, os.path.dirname(filename)) tile_layout = tiles(cropped_size, tile_size) #swapsie = tile_layout[0] #tile_layout[0] = tile_layout[1] #tile_layout[1] = swapsie tile_sizes = [ int(math.ceil(x)) for x in [ cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1] ] ] print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name) for t_y in range(tile_layout[1]-1): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) print(tile) if tile in exclude: 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] - 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) print (src_corner) bounds = create_tile(img, "%s/%s" % (directory, outfile), src_corner, src_size, quality) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile) 
import os, math, re, sys, subprocess, logging from optparse import OptionParser from osgeo import gdal import os.path """ ########################################################################## Variable assignment """ directory = 'C:\\Temp\\' #'--dir', dest='directory', help='Where to create jpeg tiles' border = 0 #'--crop', default=0, dest='border', type='int', help='Crop border' name = 'BigBoyMap' #'--name', dest='name', help='KML folder name for output' order = 20 #'--draw-order', dest='order', type='int', default=20, help='KML draw order' tile_size = 1024 #'--tile-size', dest='tile_size', default=1024, type='int', help='Max tile size [1024]' quality = 90 #'--quality', dest='quality', default=75, type='int', help='JPEG quality [75]' verbose = True #'--verbose', dest='verbose', action='store_true', help='Verbose output') extension = '.tif' compressOptions = 'COMPRESS=LZW|PREDICTOR=2|NUM_THREADS=ALL_CPUS|BIGTIFF=IF_NEEDED|TILED=YES' """ ########################################################################## Prepping of image for the tiling part """ dest = directory + name + '.kml' upResName = name + '_ur' wgsName = upResName + '_rp' finalName = wgsName + '_3b' source = directory + finalName + extension #Get the pixel size of the input raster for bumping up the res inputRaster = QgsRasterLayer(directory + name + extension) pixelSizeX = inputRaster.rasterUnitsPerPixelX() pixelSizeY = inputRaster.rasterUnitsPerPixelY() #Bump up resolution before the reproject processing.run("gdal:warpreproject", {'INPUT':directory + name + extension, 'SOURCE_CRS':None,'TARGET_CRS':None,'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-r cubic -tr ' + str(pixelSizeX / 1.5) + ' ' + str(pixelSizeY / 1.5), 'OUTPUT':directory + upResName + '.tif'}) #We need to reproject into WGS84 processing.run("gdal:warpreproject", {'INPUT':directory + upResName + extension, 'SOURCE_CRS':None,'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'RESAMPLING':0,'NODATA':255, 'TARGET_RESOLUTION':None,'OPTIONS':compressOptions,'DATA_TYPE':0,'TARGET_EXTENT':None,'TARGET_EXTENT_CRS':None,'MULTITHREADING':True, 'EXTRA':'-srcnodata 255 -dstnodata 255 -nosrcalpha -r nearest', 'OUTPUT':directory + wgsName + '.tif'}) #Remove the alpha band processing.run("gdal:translate", {'INPUT':directory + wgsName + '.tif','TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'),'NODATA':None, 'COPY_SUBDATASETS':False,'OPTIONS':compressOptions, 'EXTRA':'-a_nodata none -b 1 -b 2 -b 3','DATA_TYPE':0,'OUTPUT':source}) """ ########################################################################## Begin the function creation and the rest of the process """ 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 ] 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] ] #Create a jpeg of the given area and return the bounds. def create_tile(source, filename, offset, size, quality=75): mem_drv = gdal.GetDriverByName('MEM') mem_ds = mem_drv.Create('', size[0], size[1], source.RasterCount) bands = list(range(1, source.RasterCount+1)) data = source.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) 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(xxx_todo_changeme): (x, y) = xxx_todo_changeme 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] ]) result = { 'north': nw[1], 'east': se[0], 'south': se[1], 'west': nw[0], } jpeg_ds = None mem_ds = None return result #Create a kml file and associated images for the given georeferenced image def create_kml(source, filename, directory, tile_size=1024, border=0, name=None, order=20, exclude=[], quality=75): img = gdal.Open(source) 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 path = os.path.relpath(directory, os.path.dirname(filename)) tile_layout = tiles(cropped_size, tile_size) #swapsie = tile_layout[0] #tile_layout[0] = tile_layout[1] #tile_layout[1] = swapsie tile_sizes = [ int(math.ceil(x)) for x in [ cropped_size[0] / tile_layout[0], cropped_size[1] / tile_layout[1] ] ] print(tile_sizes) logging.debug('Using tile layout %s -> %s' % (tile_layout, tile_sizes)) bob = open(filename, 'w') bob.write("""<?xml version="1.0" encoding="UTF-8"?> <kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom"> <Folder> <name>%s</name> """ % name) for t_y in range(tile_layout[1]-1): for t_x in range(tile_layout[0]): tile = "%d,%d" % (t_y, t_x) logging.debug(tile) print(tile) if tile in exclude: 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] - 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) print (src_corner) bounds = create_tile(img, "%s/%s" % (directory, outfile), src_corner, src_size, quality) bob.write(""" <GroundOverlay> <name>%s</name> <color>ffffffff</color> <drawOrder>%d</drawOrder> <Icon> <href>%s/%s</href> <viewBoundScale>0.75</viewBoundScale> </Icon> <LatLonBox> """ % (outfile, order, path, outfile)) bob.write(""" <north>%(north)s</north> <south>%(south)s</south> <east>%(east)s</east> <west>%(west)s</west> <rotation>0</rotation> """ % bounds) bob.write(""" </LatLonBox> </GroundOverlay> """); bob.write(""" </Folder> </kml> """) bob.close() img = None """ ########################################################################## Starting up the functions """ if verbose: logging.basicConfig(level=logging.DEBUG) # set the default folder for jpegs if not directory: directory = "%s.files" % os.path.splitext(dest)[0] if not os.path.exists(directory): os.mkdir(directory) logging.info('Writing jpegs to %s' % directory) # load the exclude file exclude_file = source + ".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, directory, tile_size, border, finalName, order, exclude, quality) """ ########################################################################## Converting the kml to kmz """ def htc(m): return chr(int(m.group(1),16)) def urldecode(url): rex=re.compile('%([0-9a-hA-H][0-9a-hA-H])',re.M) return rex.sub(htc,url) #In and out paths kmlPath = dest outfile = kmlPath.replace('kml', 'kmz') if verbose: logging.basicConfig(level=logging.DEBUG) # check for KML if len(kmlPath)<1: parser.error('Path to KML file is required') if not os.path.exists(kmlPath): parser.error('Unable to find KML: %s' % kmlPath) # create the output zip file import zipfile zip = zipfile.ZipFile(outfile, 'w', zipfile.ZIP_DEFLATED) # read the source xml from xml.dom.minidom import parse kml = parse(kmlPath) nodes = kml.getElementsByTagName('href') base = os.path.dirname(directory) for node in nodes: href = node.firstChild img = urldecode(href.nodeValue).replace('file:///', '') if not os.path.exists(img): img = base + '/' + img if not os.path.exists(img): parser.error('Unable to find image: %s' % img) # add the image filename = 'files/%s' % os.path.basename(img) logging.debug("Storing %s as %s" % (img, filename)) zip.write(img, filename, zipfile.ZIP_STORED) # modify the xml to point to the correct image href.nodeValue = filename logging.debug("Storing KML as doc.kml") zip.writestr('doc.kml', kml.toxml("UTF-8")); zip.close() logging.info("Finished") print('Yeah all done, have a look at ' + outfile) 
added 10058 characters in body
Source Link
Buff Fox
  • 305
  • 1
  • 14
Loading
Source Link
Buff Fox
  • 305
  • 1
  • 14
Loading