#!/usr/bin/env python import sys import gdal import gdalconst import os import Numeric class chunk: # Create a chunk given its x, y position in the source image. # The returned image will extend to the lowest rightest pixel in both # the source and chunk image. def __init__(self, source_image, source_x, source_y): #chunk_size_X = chunk_size_Y = 64000.0 # 1000 pixels/side - GDAL barfs #chunk_size_X = chunk_size_Y = 128000.0 # 2000 pixels/side - GDAL barfs #chunk_size_X = chunk_size_Y = 256000.0 # 4000 pixels/side - GDAL barfs chunk_size_X = chunk_size_Y = 320000.0 # 5000 pixels/side chunks_path = './chunked-new' if not os.path.exists(chunks_path): os.system('mkdir -p %s' % (chunks_path)) print 'chunk(%d, %d)' % (source_x, source_y) source_geotransform = source_image.GetGeoTransform() source_projection = source_image.GetProjection() # We assume that source and chunk pixel sizes are the same. pixel_size_X = source_geotransform[1] pixel_size_Y = source_geotransform[5] print 'pixel size:', pixel_size_X, pixel_size_Y # round/floor/int? source_X = int(source_geotransform[0] + pixel_size_X * source_x) source_Y = int(source_geotransform[3] + pixel_size_Y * source_y) print 'source_location:', source_X, source_Y chunk_left_X = int(divmod(source_X, chunk_size_X)[0]) * chunk_size_X chunk_upper_Y = int(divmod(source_Y, -chunk_size_Y)[0]) * -chunk_size_Y print 'chunk image UL corner = (%f, %f)' % (chunk_left_X, chunk_upper_Y) self.chunk_basename = '%s/chunk_%d_%d' % (chunks_path, chunk_left_X / chunk_size_X, chunk_upper_Y / chunk_size_Y) self.chunk_filename = self.chunk_basename + '.tif' print 'chunk_filename=%s' % (self.chunk_filename) # If this is a new chunk, we need to create it. if not os.path.exists(self.chunk_filename): chunk_right_X = chunk_left_X + chunk_size_X chunk_lower_Y = chunk_upper_Y - chunk_size_Y chunk_size_x = int(chunk_size_X / pixel_size_X) chunk_size_y = int(chunk_size_Y / -pixel_size_Y) print 'chunk_size = (%d, %d)' % (chunk_size_x, chunk_size_y) self.image = gdal.GetDriverByName('GTiff').Create( self.chunk_filename, chunk_size_x, chunk_size_y, #5, gdal.GDT_Byte, # R,G,B,mask,sectional 4, gdal.GDT_Byte, [ 'TILED=YES' ] ) self.image.SetProjection(source_projection) self.geotransform = (chunk_left_X, pixel_size_X, 0.0, chunk_upper_Y, 0.0, pixel_size_Y) self.image.SetGeoTransform(self.geotransform) # Turn on (white) R, G and B bands. on_block = Numeric.ones((chunk_size_y, chunk_size_x)) * 256 self.image.GetRasterBand(1).WriteArray(on_block, 0, 0) self.image.GetRasterBand(2).WriteArray(on_block, 0, 0) self.image.GetRasterBand(3).WriteArray(on_block, 0, 0) # Zero the mask. off_block = Numeric.zeros((chunk_size_y, chunk_size_x)) self.image.GetRasterBand(4).WriteArray(off_block) else: # Just read an existing chunk. self.image = gdal.Open(self.chunk_filename, gdalconst.GA_Update) self.geotransform = self.image.GetGeoTransform() print 'geotransform:', self.geotransform self.ul_x = int((source_X - chunk_left_X) / pixel_size_X) self.ul_y = int((source_Y - chunk_upper_Y) / pixel_size_Y) print 'chunk\'s ul pixel:', self.ul_x, self.ul_y self.width = min(self.image.RasterXSize - self.ul_x, source_image.RasterXSize - source_x) self.height = min(self.image.RasterYSize - self.ul_y, source_image.RasterYSize - source_y) print 'size=%d, %d' % (self.width, self.height) self.band_data = [ self.image.GetRasterBand(1).ReadAsArray(self.ul_x, self.ul_y, self.width, self.height), self.image.GetRasterBand(2).ReadAsArray(self.ul_x, self.ul_y, self.width, self.height), self.image.GetRasterBand(3).ReadAsArray(self.ul_x, self.ul_y, self.width, self.height), self.image.GetRasterBand(4).ReadAsArray(self.ul_x, self.ul_y, self.width, self.height), ] def __del__(self): print 'closing' for band in range(len(self.band_data)): self.image.GetRasterBand(band + 1).WriteArray(self.band_data[band], self.ul_x, self.ul_y) self.image = None return(0) command = 'convert -scale 100%% "%s" "%s.jpg"' % (self.chunk_filename, self.chunk_basename) print command os.system(command) command = 'convert -scale 1%% "%s" "%s_tn.jpg"' % (self.chunk_filename, self.chunk_basename) print command os.system(command) command = 'convert -scale 10%% "%s" "%s_tn.png"' % (self.chunk_filename, self.chunk_basename) print command os.system(command) def warped_image_to_chunks(warped_filename): print "warped_filename=%s" % (warped_filename) warped_image = gdal.Open(warped_filename, gdalconst.GA_ReadOnly) # The upper left corner (pixel) of the area we're going to digest next. AOI_x = AOI_y = 0 if 0: # Force this to be what we want. (X, Y) = (28 * chunk_size_X, -7 * chunk_size_Y) (a, b, c, d, e, f) = warped_image.GetGeoTransform() AOI_x = ((c * (Y - d) - f * (X - a)) / (c * e - b * f)) AOI_y = ((b * (Y - d) - e * (X - a)) / (b * f - c * e)) while 1: print 'AOI: (%d, %d)' % (AOI_x, AOI_y) # Get a chunk for this area. # Clean out the old chunk before we start a new one. # If somehow we'd manage to open the same one, Bad Things could happen. chunk_AOI = None chunk_AOI = chunk(warped_image, AOI_x, AOI_y) if chunk_AOI.height == 0 or chunk_AOI.width == 0: raise 'zero dimension (%f, %f)' % (chunk_AOI.width, chunk_AOI_height) else: # Get the same region from the warped image. warped_band_data = [ warped_image.GetRasterBand(1).ReadAsArray(AOI_x, AOI_y, chunk_AOI.width, chunk_AOI.height), warped_image.GetRasterBand(2).ReadAsArray(AOI_x, AOI_y, chunk_AOI.width, chunk_AOI.height), warped_image.GetRasterBand(3).ReadAsArray(AOI_x, AOI_y, chunk_AOI.width, chunk_AOI.height), warped_image.GetRasterBand(4).ReadAsArray(AOI_x, AOI_y, chunk_AOI.width, chunk_AOI.height), ] # Update the chunk with the warped image. chunk_mask = Numeric.array(chunk_AOI.band_data[3], 'i') warped_mask = Numeric.array(warped_band_data[3], 'i') if 1: # - furthest from edge algorithm - # Choose points with the highest "distance to edge" (mask) # value and use only those. There is no blending/feathering. # These are the points we want to change. #update_mask = Numeric.greater(warped_mask, chunk_mask) # Apparently I need a bit of inertia to get rid of "streaks." update_mask = Numeric.greater(warped_mask - chunk_mask, 10) # For these "higher quality" points, copy them from the warped image. for band in range(4): Numeric.putmask(chunk_AOI.band_data[band], update_mask, warped_band_data[band]) #Numeric.putmask(chunk_AOI.band_data[band], update_mask, update_mask) elif 1: # - simple blending algorithm # Weight all values by their "distance to edge" (mask) values, # then sum and divide by weighting. # If more than two images overlap, subsequent images will seem # to be weighted more because they are averaged with the average # of the proceeding images. # Total of the weights or 1 if neither mask is set. total_weighting = Numeric.maximum(warped_mask + chunk_mask, 1) for band in range(3): chunk_AOI.band_data[band] = ( chunk_AOI.band_data[band] * chunk_mask + warped_band_data[band] * warped_mask ) / total_weighting # Set new mask to the best of the masks used. chunk_AOI.band_data[3] = Numeric.maximum(warped_mask, chunk_mask) # Set up for the next chunk. AOI_x += chunk_AOI.width # Are we right of the right edge? if AOI_x >= warped_image.RasterXSize: AOI_x = 0 print 'skipping down %d' % (chunk_AOI.height) AOI_y += chunk_AOI.height # Are we on or below the bottom edge? if AOI_y >= warped_image.RasterYSize: print 'breaking on ', AOI_y, warped_image.RasterXSize break # break for filename in sys.argv[1:]: warped_image_to_chunks(filename)