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 a geotiff image with pixel intensities ranging between certain values. Lets say bewteen -2 and 2 (the image is normalized). I have a certain threshold which interests me, lets sat 1. I would like to create a NEW image that makes all the pixels that where previously larger than 1 to now have the maximum value, 2, or just some constant value. And all other pixels to be transparent. Im very new to this so I hope my question is understandable.

greetings and thank you in advanced.

share|improve this question
What software are you using? Is this an ESRI question with Spatial Analyst, or an Open Source, QGIS, question? –  Get Spatial Jun 14 '12 at 4:52
What is ESRI? I'm sorry. I just recently started working on these things. At the moment I just play around on images with python and then visualize my results using ENVI. What program would you recommend, by the way, for visualization and such? I don't think im loving ENVI. –  JEquihua Jun 14 '12 at 18:57
ESRI is a company that makes a commercial GIS package called ArcGIS. It is a good product, but rather expensive. An open source option for visualizing raster data is Quantum GIS. –  Get Spatial Jun 14 '12 at 20:05
You'll find a ton of information in those slides: gis.usu.edu/~chrisg/python/2009 –  WAF Aug 28 '13 at 12:43
add comment

1 Answer

up vote 5 down vote accepted

If you're using Python I'd recommend using the GDAL library, which has it's own Python bindings. Assuming you've got both GDAl (see this GIS StackExchange question for details on how to install on windows) and numpy installed, your code could look something like:

from osgeo import gdal
import numpy as np

#Open our original data as read only
dataset = gdal.Open("file_path.tif", gdal.GA_ReadOnly)

#Note that unlike the rest of Python, Raster Bands in GDAL are numbered
#beginning with 1.
#I suspect this is to conform to the landsat band naming convention
band = dataset.GetRasterBand(1)

#Read in the data from the band to a numpy array
data = band.ReadAsArray()
data = data.astype(numpy.float)

#Use numpy, scipy, and whatever Python to make some output data
#That for ease of use should be an array of the same size and dimensions
#as the input data.
out_data = np.where(abs(data) > 1., 2., -9999.)
#Note -9999 is a convenience value for null - there's no number for
#transparent values - it's just how you visualise the data in the viewer

#And now we start preparing our output
driver = gdal.GetDriverByName("GTiff")
metadata = driver.GetMetadata()

#Create an output raster the same size as the input
out = driver.Create("out_file.tif",
                    1, #Number of bands to create in the output

#Copy across projection and transform details for the output

#Get the band to write to
out_band = out.GetRasterBand(1)

#And write our processed data
share|improve this answer
This looks amazing, I suspect just what I need. Let me try it out and get back to you. By the way, if there is no number for transparent values I will have to indicate my visualization (e.g. ENVI) that I want X intensity values to be transparent? @om_henners –  JEquihua Jun 14 '12 at 18:53
That's right @JEquihua - most visualisation tools have the option of setting values transparent. –  om_henners Jun 16 '12 at 1:18
I noticed that you're assigning "metadata", but don't seem to be doing anything with it. What does "metadata" provide, and how would you use it? –  monkut Jun 18 '12 at 4:41
add comment

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.