import sys import numpy import commands import optparse from osgeo import gdal, osr class Quadrangle: def __init__(self, filename): self.filename = filename self.dem = gdal.Open(self.filename) assert self.dem.RasterCount == 1 # my kingdom for GetUnitType in osgeo.gdal status, output = commands.getstatusoutput('gdalinfo "%s" | grep "Unit Type:"' % self.filename) assert status == 0 self.units = output.split(': ')[1] assert self.units in ('ft', 'm') self.srs = osr.SpatialReference() self.srs.ImportFromWkt(self.dem.GetProjectionRef()) self.band = self.dem.GetRasterBand(1) self.nodata = self.band.GetNoDataValue() tx = self.dem.GetGeoTransform() # zero rotations assert tx[2] == 0 and tx[4] == 0 self.xoffset, self.yoffset = tx[0], tx[3] self.xscale, self.yscale = tx[1], tx[5] self.width, self.height = self.dem.RasterXSize, self.dem.RasterYSize def pixelPoint(self, x, y): """ Return (x, y) point in projected space for a given (x, y) image pixel. """ return (x * self.xscale) + self.xoffset, (y * self.yscale) + self.yoffset def pointPixel(self, x, y): """ Return (x, y) image pixel for a given (x, y) point in projected space. """ return (x - self.xoffset) / self.xscale, (y - self.yoffset) / self.yscale def pointIncluded(self, x, y): """ True if given point in projected space exists as an image pixel. """ x, y = self.pointPixel(x, y) return 0 <= x and x < self.width and 0 <= y and y < self.height def isCompatibleWith(self, other): """ True if other quadrangle matches local projection and scale. """ return (self.xscale, self.yscale) == (other.xscale, other.yscale) \ and self.srs.ExportToWkt() == other.srs.ExportToWkt() def getArray(self): """ Get a new numpy array of raster data as float32. """ data = self.band.ReadRaster(0, 0, self.width, self.height, buf_type=gdal.GDT_Float32) array = numpy.fromstring(data, dtype=numpy.float32).reshape(self.height, self.width) if self.units == 'ft': # optionally convert to meters defined = array != self.nodata array[defined] = array[defined] / 3.2808399 return array def adjustSlices(aTopleft, aBottomright, bTopleft, bBottomright, aWidth, aHeight, bWidth, bHeight): """ """ aTopleft_ = list(aTopleft) aBottomright_ = list(aBottomright) bTopleft_ = list(bTopleft) bBottomright_ = list(bBottomright) if aTopleft_[0] < 0: diff = -aTopleft_[0] aTopleft_[0] += diff bTopleft_[0] += diff if aTopleft_[1] < 0: diff = -aTopleft_[1] aTopleft_[1] += diff bTopleft_[1] += diff if bTopleft_[0] < 0: diff = -bTopleft_[0] aTopleft_[0] += diff bTopleft_[0] += diff if bTopleft_[1] < 0: diff = -bTopleft_[1] aTopleft_[1] += diff bTopleft_[1] += diff if aBottomright_[0] > aWidth: diff = aBottomright_[0] - aWidth aBottomright_[0] -= diff bBottomright_[0] -= diff if aBottomright_[1] > aHeight: diff = aBottomright_[1] - aHeight aBottomright_[1] -= diff bBottomright_[1] -= diff if bBottomright_[0] > bWidth: diff = bBottomright_[0] - bWidth aBottomright_[0] -= diff bBottomright_[0] -= diff if bBottomright_[1] > bHeight: diff = bBottomright_[1] - bHeight aBottomright_[1] -= diff bBottomright_[1] -= diff return aTopleft_, aBottomright_, bTopleft_, bBottomright_ if __name__ == '__main__': parser = optparse.OptionParser('usage: ...') parser.set_defaults(buffer=100, verbose=False) parser.add_option('-v', '--verbose', dest='verbose', action='store_true', help='Verbose output') parser.add_option('-b', '--buffer', dest='buffer', type='int', help='Desired buffer in pixels') parser.add_option('-c', '--center', dest='center', help='Center .dem file') parser.add_option('-o', '--output', dest='output', help='Output .dem file') (options, args) = parser.parse_args() buffer = options.buffer center = Quadrangle(options.center) output = options.output neighbors = [neighbor for neighbor in map(Quadrangle, args) if neighbor.isCompatibleWith(center)] centerarray = numpy.empty((center.height + buffer*2, center.width + buffer*2), dtype=numpy.float32) centerarray[:, :] = center.nodata centerarray[buffer:-buffer, buffer:-buffer] = center.getArray() # clockwise from top-left pixels = (('up left', -1, -1, buffer, buffer), ('up', center.width/2, -1, center.width/2, buffer), ('up right', center.width+1, -1, buffer, buffer), ('right', center.width+1, center.height/2, buffer, center.height/2), ('down right', center.width+1, center.height+1, buffer, buffer), ('down', center.width/2, center.height+1, center.width/2, buffer), ('down left', -1, center.height+1, buffer, buffer), ('left', -1, center.height/2, buffer, center.height/2)) # pixels = ((center.width+1, -1, buffer, buffer), ) for (direction, x, y, w, h) in pixels: point = center.pixelPoint(x, y) for neighbor in neighbors: if neighbor.pointIncluded(*point): if options.verbose: print >> sys.stderr, neighbor.filename, direction, 'of', center.filename ctopleft = [x - w, y - h] cbottomright = [x + w, y + h] ntopleft = map(round, neighbor.pointPixel(*center.pixelPoint(*ctopleft))) nbottomright = map(round, neighbor.pointPixel(*center.pixelPoint(*cbottomright))) ctopleft = [ctopleft[0] + buffer, ctopleft[1] + buffer] cbottomright = [cbottomright[0] + buffer, cbottomright[1] + buffer] # slices must only cover existing data ctopleft, cbottomright, ntopleft, nbottomright \ = adjustSlices(ctopleft, cbottomright, ntopleft, nbottomright, centerarray.shape[1], centerarray.shape[0], neighbor.width, neighbor.height) # reference just the portions of center and neighbor arrays we care about dst = centerarray[ctopleft[1]:cbottomright[1], ctopleft[0]:cbottomright[0]] src = neighbor.getArray()[ntopleft[1]:nbottomright[1], ntopleft[0]:nbottomright[0]] # assign data where none exists dst[dst == center.nodata] = src[dst == center.nodata] driver = gdal.GetDriverByName('GTiff') out = driver.Create(output, centerarray.shape[1], centerarray.shape[0], 1, gdal.GDT_Int16) outtx = list(center.dem.GetGeoTransform()) outtx[0], outtx[3] = center.pixelPoint(-buffer, -buffer) out.SetGeoTransform(outtx) out.SetProjection(center.dem.GetProjection()) outband = out.GetRasterBand(1) outband.SetNoDataValue(center.band.GetNoDataValue()) outdata = centerarray.astype(numpy.int16).tostring() outband.WriteRaster(0, 0, out.RasterXSize, out.RasterYSize, outdata, buf_type=gdal.GDT_Int16) sys.exit(0)