Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have small program in Python that uses gdal and numpy to perform a logarithmic scale of some data into the RGB 0-255 range.

Now this input data has a nodata value set to mask the valid values, and I would like to have the output to retain that mask.

Obviously I can't use nodata for my output as I am ussing the full byte range for colour values.

From what I can tell from the GDAL API I need to create a mask band in the output file.

I added the code to do this but suddenly the program started running really slow and the output file size became huge.

What it the best way to do this?

Here is my code:

#!/usr/bin/env python

import gdal
from gdalconst import *
import numpy
import sys

infile = sys.argv[1]
outfile = sys.argv[2]
format = 'GTiff'
type = GDT_Byte

indataset = gdal.Open(infile, GA_ReadOnly)

gdal.SetConfigOption('GDAL_TIFF_INTERNAL_MASK', 'YES')
out_driver = gdal.GetDriverByName(format)
outdataset = out_driver.Create(outfile, indataset.RasterXSize, indataset.RasterYSize, indataset.RasterCount, type)
outdataset.CreateMaskBand(gdal.GMF_PER_DATASET)

gt = indataset.GetGeoTransform()
outdataset.SetGeoTransform(gt)

prj = indataset.GetProjectionRef()
outdataset.SetProjection(prj)

for iBand in range(1, indataset.RasterCount + 1):
    inband = indataset.GetRasterBand(iBand)
    inmaskband = inband.GetMaskBand()

    outband = outdataset.GetRasterBand(iBand)
    outmaskband = outband.GetMaskBand()

    for i in range(inband.YSize - 1, -1, -1):
        scanline = inband.ReadAsArray(0, i, inband.XSize, 1, inband.XSize, 1)
        scanline = numpy.multiply(numpy.divide(numpy.log1p(numpy.minimum(10000, numpy.maximum(0, scanline))), numpy.log1p(10000)), 255)
        outband.WriteArray(scanline, 0, i)

        scanline = inmaskband.ReadAsArray(0, i, inband.XSize, 1, inband.XSize, 1)
        outmaskband.WriteArray(scanline, 0, i)
share|improve this question
    
Why are you creating masks from all bands? Did you tried this gis.stackexchange.com/questions/27426/… –  rashad Nov 2 '12 at 5:14
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.