#!/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 import math 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 = warped_y_pixel_size = 1.0 / pixels_per_degree 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), ((x_min + x_max) / 2, y_min), ((x_min + x_max) / 2, y_max), (x_min, (y_min + y_max) / 2), (x_max, (y_min + y_max) / 2), ] 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) * pixels_per_degree warped_ysize = (targ_ymax - targ_ymin) * pixels_per_degree 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 60%% "%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 = degrees_per_chunk chunks_path = './data/square_chunked' 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 source_X = source_geotransform[0] + pixel_size_X * source_x source_Y = source_geotransform[3] + pixel_size_Y * source_y print 'source_location:', source_X, source_Y # Chunks are named by the left upper corners. chunk_left_X = int(math.floor(source_X)) chunk_upper_Y = int(math.ceil(source_Y)) print 'chunk image left uppercorner = (%d, %d)' % (chunk_left_X, chunk_upper_Y) self.chunk_basename = '%s/chunk_%d_%d' % (chunks_path, chunk_left_X, chunk_upper_Y) self.chunk_filename = self.chunk_basename + '.tif' print 'chunk_filename=%s' % (self.chunk_filename) self.update = True # If this is a new chunk, we need to create it. if not os.path.exists(self.chunk_filename): self.new_chunk = True chunk_size_x = chunk_size_y = pixels_per_degree * degrees_per_chunk 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: self.new_chunk = False # Just read an existing chunk. self.image = gdal.Open(self.chunk_filename, gdalconst.GA_Update) self.geotransform = self.image.GetGeoTransform() print 'geotransform:', self.geotransform print 'source_X=%f, source_Y=%f' % (source_X, source_Y) print 'chunk_left_X=%f, chunk_upper_Y=%f' % (chunk_left_X, chunk_upper_Y) self.ul_x = int((source_X - chunk_left_X) * pixels_per_degree) self.ul_y = int(-(source_Y - chunk_upper_Y) * pixels_per_degree) #self.ul_x = ((source_X - chunk_left_X) * pixels_per_degree) #self.ul_y = (-(source_Y - chunk_upper_Y) * pixels_per_degree) #self.ul_y = (chunk_upper_Y - source_Y) * pixels_per_degree 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 disable_update(self): print 'disable_update' self.update = False def __del__(self): print 'closing' # Are we updating this chunk? if self.update: print 'writing chunk' for band in range(len(self.band_data)): self.image.GetRasterBand(band + 1).WriteArray(self.band_data[band], self.ul_x, self.ul_y) elif self.new_chunk: print 'removing newly-created chunk file' self.image = None os.unlink(self.chunk_filename) self.image = None return(0) command = 'convert -scale 100%% "%s" "%s.jp2"' % (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) print 'chunk size: (%d, %d)' % (chunk_AOI.width, chunk_AOI.height) if chunk_AOI.height == 0 or chunk_AOI.width == 0: #continue 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: # - further 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) if update_mask: print 'updating' # 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) else: chunk_AOI.disable_update() else: # - 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 at 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/square_warped/' # Load info about the center of our universe. reference_projection = 'GEOGCS["unnamed",DATUM["unknown",SPHEROID["unnamed",6378137,298.2572235629972]],PRIMEM["Greenwich",0],UNIT[,0.0174532925199433]]' reference_srs = osr.SpatialReference() reference_srs.ImportFromWkt(reference_projection) 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 == 'Dallas-Ft Worth': 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) pixels_per_degree = 1350 degrees_per_chunk = 1 if 1: process_all('/var/www/aviationtoolbox/raw_data/FAA/sectionals/new/') else: reference_projection = 'GEOGCS["unnamed",DATUM["unknown",SPHEROID["unnamed",6378137,298.2572235629972]],PRIMEM["Greenwich",0],UNIT[,0.0174532925199433]]' reference_srs = osr.SpatialReference() reference_srs.ImportFromWkt(reference_projection) #warp_masked_file('data/masked/Albuquerque 73 South.tif', 'data/square_warped/Albuquerque 73 South.tif', reference_srs) #warped_image_to_chunks('data/square_warped/Albuquerque 73 South.tif') warped_image_to_chunks('data/square_warped/Chicago 68 North.tif') warped_image_to_chunks('data/square_warped/Chicago 68 South.tif')