#!/usr/bin/env python import gdal import gdalconst import sys import Numeric import os.path import osr import colorsys import Image import glob import string import os import pickle import time def mask_min(masked_band, mask, x, y): (size_y, size_x) = Numeric.shape(mask) current = masked_band.ReadAsArray(x, y, size_x, size_y) masked_band.WriteArray(Numeric.minimum(current, mask), x, y) def find_x_pos(ct_to_lonlat, geomatrix, y, target_lon, hint_x=0): (lon, lat) = transform_pixel(ct_to_lonlat, geomatrix, (hint_x, y)) x = hint_x # If we're to the right of the left edge, move left until we're # left of the edge and then add one. if lon > target_lon: while (lon > target_lon): x -= 1 (lon, lat) = transform_pixel(ct_to_lonlat, geomatrix, (x, y)) return(x + 1, lat) # If we're to the left of the left edge, move right until we're # right of the edge. else: while (lon < target_lon): x += 1 (lon, lat) = transform_pixel(ct_to_lonlat, geomatrix, (x, y)) return(x, lat) def find_y_pos(ct_to_lonlat, geomatrix, x, target_lat, hint_y=0): (lon, lat) = transform_pixel(ct_to_lonlat, geomatrix, (x, hint_y)) y = hint_y # If we're above the bottom edge, move down until we're # below the edge and then subtract one. if lat > target_lat: while (lat > target_lat): y += 1 (lon, lat) = transform_pixel(ct_to_lonlat, geomatrix, (x, y)) return(y - 1, lon) # If we're below the bottom edge, move up until we're # above the edge. else: while (lat < target_lat): y -= 1 (lon, lat) = transform_pixel(ct_to_lonlat, geomatrix, (x, y)) return(y, lon) def transform_pixel(ct, geomatrix, pixel): (pixel_x, pixel_y) = pixel X = geomatrix[0] + geomatrix[1] * pixel_x + geomatrix[2] * pixel_y Y = geomatrix[3] + geomatrix[4] * pixel_x + geomatrix[5] * pixel_y (lon, lat, height) = ct.TransformPoint(X, Y) return(lon, lat) def transform_lonlat(ct, geomatrix, lonlat): (a, b, c, d, e, f) = geomatrix (lon, lat) = lonlat (X, Y, Z) = ct.TransformPoint(lon, lat) x = int((c * (Y - d) - f * (X - a)) / (c * e - b * f)) y = int((b * (Y - d) - e * (X - a)) / (b * f - c * e)) return(x, y) def mask_sectional_file(sectional_filename, masked_filename, sectional_side, ll_corner_lonlat, reference_srs): print 'mask_sectional_file(%s, %s)' % (sectional_filename, masked_filename) if os.path.exists(masked_filename): print 'already have %s, skipping' % (masked_filename) return(True) temp_filename = masked_filename + '.tmp' masked_bands = 4 horizontal_skip_min = 850 horizontal_skip_max = 2000 top_border = 0 # Do we need this? edge_value = 255 sectional_image = gdal.Open(sectional_filename, gdalconst.GA_ReadOnly) sectional_band = sectional_image.GetRasterBand(1) sectional_color_table = sectional_band.GetRasterColorTable() sectional_geomatrix = sectional_image.GetGeoTransform() sectional_srs = osr.SpatialReference() sectional_projection = sectional_image.GetProjection() sectional_srs.ImportFromWkt(sectional_projection) srs_lonlat = sectional_srs.CloneGeogCS() sectional_ct_to_lonlat = osr.CoordinateTransformation(sectional_srs, srs_lonlat) sectional_ct_to_pixel = osr.CoordinateTransformation(srs_lonlat, sectional_srs) sectional_to_reference_ct = osr.CoordinateTransformation(sectional_srs, reference_srs) masked_image = gdal.GetDriverByName('GTiff').Create( #masked_filename, temp_filename, sectional_image.RasterXSize, sectional_image.RasterYSize, masked_bands, gdal.GDT_Byte, # R,G,B,mask [ 'TILED=YES' ] ) masked_image.SetGeoTransform(sectional_geomatrix) masked_image.SetProjection(sectional_projection) masked_band = masked_image.GetRasterBand(1) masked_band.SetDescription(sectional_band.GetDescription()) masked_band.SetMetadata(sectional_band.GetMetadata()) # masked_band.SetRasterColorTable(sectional_band.GetRasterColorTable()) # colorish_index = Numeric.ones (256) * 4 masked_band_r = masked_image.GetRasterBand(1) masked_band_g = masked_image.GetRasterBand(2) masked_band_b = masked_image.GetRasterBand(3) #masked_band_a = masked_image.GetRasterBand(4) masked_band_m = masked_image.GetRasterBand(4) # mask quality # This doesn't seem to work. #masked_band_r.SetNoDataValue(255) #masked_band_g.SetNoDataValue(255) #masked_band_b.SetNoDataValue(255) #masked_band_a.SetNoDataValue(0) #masked_band_m.SetNoDataValue(255) # Initialize the image manually. blank_lines = Numeric.ones((sectional_image.RasterYSize, sectional_image.RasterXSize)) * 255 # For testing. # masked_band_m = masked_band_r #masked_band_r.WriteArray(blank_lines, 0, 0) #masked_band_g.WriteArray(blank_lines, 0, 0) #masked_band_b.WriteArray(blank_lines, 0, 0) masked_band_m.WriteArray(blank_lines, 0, 0) # Create the palette. lookup = [ Numeric.arrayrange(257), Numeric.arrayrange(257), Numeric.arrayrange(257), Numeric.ones(257)*255 ] for i in range (min (256, sectional_color_table.GetCount ())): #print i, (r, g, b, x) = sectional_color_table.GetColorEntry (i) (h, l, s) = colorsys.rgb_to_hls(r, g, b) #print '%0.3f, %0.3f, %03d' % (h, l, s), if s > -0.04 and l > 128: colorish_index[i] = 0 else: colorish_index[i] = 1 # If it's completely white, pull a dash of blue. #if r + g + b == 3 * 255: # b = 254 lookup[0][i] = r lookup[1][i] = g lookup[2][i] = b # I think this is a mask for the edge. colorish_index[edge_value] = -5 # Where is the lower left corner? ll_corner_pixel = transform_lonlat(sectional_ct_to_pixel, sectional_geomatrix, ll_corner_lonlat) # We'll use these a lot. (ll_corner_x, ll_corner_y) = ll_corner_pixel (ll_corner_lon, ll_corner_lat) = ll_corner_lonlat (left_edge_x, lat) = find_x_pos(sectional_ct_to_lonlat, sectional_geomatrix, 0, ll_corner_lon, hint_x=ll_corner_x) x_min = x_max = sectional_image.RasterXSize y_min = y_max = 0 # If it's the North side, skip down until we're in colorful lines. if sectional_side == 'North': # Find the top line. inmap_vertical_count = 0 for iY in range(top_border, sectional_image.RasterYSize, 1): sectional_data = sectional_band.ReadAsArray( ll_corner_x, iY, sectional_image.RasterXSize - ll_corner_x, 1 ) sectional_data_indexed = Numeric.take(colorish_index, sectional_data)[0] colorval = sum(sectional_data_indexed) * 100 / len(sectional_data_indexed) print colorval # print colorval #if colorval >= 80: if colorval >= 58: inmap_vertical_count += 1 if inmap_vertical_count > 3: y_min = iY # - 3 break else: inmap_vertical_count = 0 # For South sides, we need to feel for the top "line" - it might not be # horizontal. (See Great Falls.) else: # y_min = top_border print 'Feeling for the top edge.' mask = Numeric.fromfunction(lambda y, x: y, (255, 1)) # Detect the top edge. # hardcoded vertical extent for iX in range(left_edge_x, 1000): for iY in range(0, sectional_image.RasterYSize): sectional_band_data = sectional_band.ReadAsArray(iX, iY, 1, 1) if not sectional_band_data[0][0] == edge_value: break #print iY y_min = min(y_min, iY) # Mask out above top edge. masked_band_m.WriteArray(Numeric.zeros((iY, 1)), iX, 0) # Feather mask below bottom edge. mask_min(masked_band_m, mask, iX, iY) # Mask out everything above top line as 0. masked_band_m.WriteArray(Numeric.zeros((y_min, sectional_image.RasterXSize)), 0, 0) # Feather mask below the top line. mask = Numeric.fromfunction(lambda y, x: y, (255, sectional_image.RasterXSize)) # It's the first one, so don't bother checking what's there. masked_band_m.WriteArray(mask, 0, y_min) # Start at the x position of the LL corner just to get close. # Use the one from top edge search. # left_edge_x = ll_corner_x mask = Numeric.fromfunction(lambda y, x: x, (1, 255)) # Find the left edge. for iY in range(y_min, sectional_image.RasterYSize, 1): # Use the previous left_edge_x as a starting point. (left_edge_x, lat) = find_x_pos(sectional_ct_to_lonlat, sectional_geomatrix, iY, ll_corner_lon, hint_x=left_edge_x) # If we're below the LL corner, stop. if lat < ll_corner_lat: break x_min = min(left_edge_x, x_min) # Mask out to the left of the left edge. masked_band_m.WriteArray(Numeric.zeros((1, left_edge_x)), 0, iY) # Feather mask to the right of the left_edge. mask_min(masked_band_m, mask, left_edge_x, iY) # If we just hit the bottom of the image, feather mask above bottom edge. if lat >= ll_corner_lat: #print 'Feeling for the bottom edge. (x_min=%d)' % (x_min) mask = Numeric.fromfunction(lambda y, x: 255 - y, (255, 1)) # Detect the bottom edge. for iX in range(x_min, sectional_image.RasterXSize): for iY in range(sectional_image.RasterYSize - 1, 255, -1): sectional_band_data = sectional_band.ReadAsArray(iX, iY, 1, 1) if not sectional_band_data[0][0] == edge_value: break y_max = max(y_max, iY) # Mask out below bottom edge. masked_band_m.WriteArray(Numeric.zeros((sectional_image.RasterYSize - iY, 1)), iX, iY) # Feather mask above bottom edge. mask_min(masked_band_m, mask, iX, iY - 255) # Otherwise we are below the LL corner now. else: # Mask out left and below the LL corner. masked_band_m.WriteArray(Numeric.zeros((sectional_image.RasterYSize - ll_corner_y, ll_corner_x)), 0, ll_corner_y) # Follow the bottom edge to the right. bottom_edge_y = ll_corner_y for iX in range(ll_corner_x, sectional_image.RasterXSize): (bottom_edge_y, lon) = find_y_pos(sectional_ct_to_lonlat, sectional_geomatrix, iX, ll_corner_lat, hint_y=bottom_edge_y) y_max = max(bottom_edge_y, y_max) # Mask out below. masked_band_m.WriteArray(Numeric.zeros((sectional_image.RasterYSize - bottom_edge_y, 1)), iX, bottom_edge_y) # Feather mask above. mask = Numeric.fromfunction(lambda y, x: 255 - y, (255, 1)) mask_min(masked_band_m, mask, iX, bottom_edge_y - 255) # Detect right edge. # I'm tempted to do this when copying all of the image data, but for now # it's here so I can tweak it more easily. print 'detecting right edge...' mask = Numeric.fromfunction(lambda y, x: 255 - x, (1, 255)) for iY in range(y_min, y_max): for iX in range(sectional_image.RasterXSize - 1, 0, -1): # Read one pixel at a time. # Is this horrible? sectional_band_data = sectional_band.ReadAsArray(iX, iY, 1, 1) if not sectional_band_data[0][0] == edge_value: break #print 'edge detect' #print iX # If it's reasonable... if iX > 255: # Mask out to the right of the right edge. masked_band_m.WriteArray(Numeric.zeros((1, sectional_image.RasterXSize - iX)), iX, iY) # Feather mask to the left of the right edge. mask_min(masked_band_m, mask, iX - 255, iY) # Otherwise, obliterate the whole band. else: masked_band_m.WriteArray(Numeric.zeros((1, sectional_image.RasterXSize)), 0, iY) # Copy over the sectional image, making completely masked (0) pixels white. # I suppose this could be done all at once instead of line by line? gdal.TermProgress(0.0) for iY in range(sectional_image.RasterYSize): #for iY in range(30): mask = masked_band_m.ReadAsArray(0, iY, sectional_image.RasterXSize, 1) sectional_band_data = sectional_band.ReadAsArray(0, iY, sectional_image.RasterXSize, 1) for band in range(3): band_lookup = lookup[band] masked_data = Numeric.array([Numeric.ones(sectional_image.RasterXSize) * 255]) Numeric.putmask(masked_data, mask, Numeric.take(band_lookup, sectional_band_data)) masked_image.GetRasterBand(band + 1).WriteArray(masked_data, 0, iY) gdal.TermProgress( (iY+1.0) / sectional_image.RasterYSize ) masked_image = None command = 'mv "%s" "%s"' % (temp_filename, masked_filename) print command os.system(command) command = 'convert -scale 60%% "%s" "%s.jpg"' % (masked_filename, masked_filename[:-4]) print command os.system(command) def warp_masked_file(masked_filename, warped_filename, reference_srs): print 'warp_masked_file(%s, %s)' % (masked_filename, warped_filename) if os.path.exists(warped_filename): print 'already have %s, skipping' % (warped_filename) return(True) temp_filename = warped_filename + '.tmp' warped_x_pixel_size = 64.0 warped_y_pixel_size = 64.0 warped_bands = 4 masked_image = gdal.Open(masked_filename, gdalconst.GA_ReadOnly) masked_geomatrix = masked_image.GetGeoTransform() masked_srs = osr.SpatialReference() masked_projection = masked_image.GetProjection() masked_srs.ImportFromWkt(masked_projection) masked_to_reference_ct = osr.CoordinateTransformation(masked_srs, reference_srs) x_min = y_min = 0 x_max = masked_image.RasterXSize y_max = masked_image.RasterYSize corners = [ (x_min, y_min), (x_min, y_max), (x_max, y_min), (x_max, y_max), ] targ_xmin = None for corner in corners: (x, y) = corner X = masked_geomatrix[0] + masked_geomatrix[1] * x + masked_geomatrix[2] * y Y = masked_geomatrix[3] + masked_geomatrix[4] * x + masked_geomatrix[5] * y (x, y, z) = masked_to_reference_ct.TransformPoint(X, Y) print corner, x,y if targ_xmin is None: targ_xmin = targ_xmax = x targ_ymin = targ_ymax = y else: targ_xmin = min(targ_xmin, x) targ_xmax = max(targ_xmax, x) targ_ymin = min(targ_ymin, y) targ_ymax = max(targ_ymax, y) #print targ_xmin, targ_xmax #print targ_ymin, targ_ymax warped_xsize = (targ_xmax - targ_xmin) / warped_x_pixel_size + 1 warped_ysize = (targ_ymax - targ_ymin) / warped_y_pixel_size + 1 print warped_xsize, warped_ysize warped_gt = (targ_xmin, warped_x_pixel_size, 0.0, targ_ymax, 0.0, - warped_y_pixel_size) warped_image = gdal.GetDriverByName('GTiff').Create( #warped_filename, temp_filename, warped_xsize, warped_ysize, warped_bands, gdal.GDT_Byte, # RGBM [ 'TILED=YES' ] ) warped_image.SetProjection( reference_srs.ExportToWkt() ) warped_image.SetGeoTransform( warped_gt ) warped_image.GetRasterBand(1).SetNoDataValue(255) warped_image.GetRasterBand(2).SetNoDataValue(255) warped_image.GetRasterBand(3).SetNoDataValue(255) warped_image.GetRasterBand(4).SetNoDataValue(255) white_lines = Numeric.ones((warped_ysize, warped_xsize)) * 255 black_lines = Numeric.zeros((warped_ysize, warped_xsize)) warped_image.GetRasterBand(1).WriteArray(white_lines,0,0) warped_image.GetRasterBand(2).WriteArray(white_lines,0,0) warped_image.GetRasterBand(3).WriteArray(white_lines,0,0) warped_image.GetRasterBand(4).WriteArray(black_lines,0,0) warped_image = None # Warp the image to match the reference. command = './warp "%s" "%s"' % (masked_filename, temp_filename) print command os.system(command) command = 'mv "%s" "%s"' % (temp_filename, warped_filename) print command os.system(command) ## Make thumbnails. #command = 'convert -scale 40%% "%s" "%s.png"' % (warped_filename, warped_filename[:-4]) #print command #os.system(command) #command = 'convert -scale 20%% "%s" "%s.jpg"' % (warped_filename, warped_filename[:-4]) #print command #os.system(command) ## Split the image into chunks. #command = './chunk.py "%s"' % (warped_filename) #print command #os.system(command) 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 = './data/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 def process_all(sectionals_path): masked_sectionals_path = './data/masked/' warped_sectionals_path = './data/warped/' # Load info about the center of our universe. #reference_filename = glob.glob('../raw_data/FAA/sectionals/new/Wichita*North.tif')[0] reference_filename = './data/reference_sectional.tif' reference_image = gdal.Open(reference_filename, gdalconst.GA_ReadOnly) reference_projection = reference_image.GetProjection() reference_srs = osr.SpatialReference() reference_srs.ImportFromWkt(reference_projection) reference_image = None time_prev = None sectionals_info_file = open('./data/sectional_corners', 'r') for sectional_info_line in sectionals_info_file.readlines(): print '*' * 80 now = time.time() if time_prev: print 'elapsed time: %d minutes' % ((now - time_prev) / 60) time_prev = now #print sectional_info_line (sectional_region, ll_corner_string, notes) = string.split(sectional_info_line, '|') if 0 and not sectional_region == 'Albuquerque': continue ll_corner_lonlat = tuple(map(float, string.split(ll_corner_string, ' '))) try: sectional_version = glob.glob('%s%s ?? North.tif' % (sectionals_path, sectional_region))[0][-12:-10] except: print 'failed to find %s sectional' % (sectional_region) continue for sectional_side in ['North', 'South']: sectional_filename = '%s%s %s %s.tif' % (sectionals_path, sectional_region, sectional_version, sectional_side) masked_filename = '%s%s %s %s.tif' % (masked_sectionals_path, sectional_region, sectional_version, sectional_side) warped_filename = '%s%s %s %s.tif' % (warped_sectionals_path, sectional_region, sectional_version, sectional_side) info_filename = 'masked/%s %s %s.info' % (sectional_region, sectional_version, sectional_side) mask_sectional_file(sectional_filename, masked_filename, sectional_side, ll_corner_lonlat, reference_srs) warp_masked_file(masked_filename, warped_filename, reference_srs) warped_image_to_chunks(warped_filename) # sys.exit(0) if 1: #process_all('/var/www/aviationtoolbox/raw_data/FAA/sectionals/new/') process_all('/var/www/aviationtoolbox/raw_data/FAA/sectionals/current/') else: sectionals_path = '/var/www/aviationtoolbox/raw_data/FAA/sectionals/new/' masked_sectionals_path = './data/masked/' warped_sectionals_path = './data/warped/' sectional_region = 'Albuquerque' sectional_version = '70' sectional_side = 'North' sectional_filename = '%s%s %s %s.tif' % (sectionals_path, sectional_region, sectional_version, sectional_side) masked_filename = '%s%s %s %s.tif' % (masked_sectionals_path, sectional_region, sectional_version, sectional_side) warped_filename = '%s%s %s %s.tif' % (warped_sectionals_path, sectional_region, sectional_version, sectional_side) info_filename = 'data/masked/%s %s %s.info' % (sectional_region, sectional_version, sectional_side) mask_sectional_file(sectional_filename, masked_filename, sectional_side, ll_corner_lonlat, reference_srs) warp_masked_file(masked_filename, warped_filename, reference_srs) warped_image_to_chunks(warped_filename)