Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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)

gt = indataset.GetGeoTransform()

prj = indataset.GetProjectionRef()

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… – rashad Nov 2 '12 at 5:14

for one particular band, here's how I create a nodata mask.

import gdal
import numpy as np

nodata = -9999#or set it based on the raster nodata value.
r = gdal.Open(rastername)
raster_arr = np.array(r.GetRasterBand(band).ReadAsArray())

nodatamask = raster_arr == nodata

#do your thing and end up with 
#a result_raster that you'd like to save

result_raster[nodatamask] = nodata
#then save.

Of course, if you're doing logarithms and you have negative values, you'd need to make a temporary nodata value.

nodatamask = raster_arr == nodata

raster_arr[nodatamask] = 100000000

#do your thing

result_raster[nodatamask] = nodata
#then save.
share|improve this answer

Your Answer


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

Not the answer you're looking for? Browse other questions tagged or ask your own question.