6

I have a GeoTIFF image with pixel intensities ranging between certain values. Let's say between -2 and 2 (the image is normalized). I have a certain threshold which interests me, let's 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.

I'm very new to this so I hope my question is understandable.

4
  • 1
    What software are you using? Is this an ESRI question with Spatial Analyst, or an Open Source, QGIS, question? Commented Jun 14, 2012 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
    Commented Jun 14, 2012 at 18:57
  • 1
    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. Commented Jun 14, 2012 at 20:05
  • You'll find a ton of information in those slides: gis.usu.edu/~chrisg/python/2009
    – WAF
    Commented Aug 28, 2013 at 12:43

1 Answer 1

6

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",
                    dataset.RasterXSize,
                    dataset.RasterYSize,
                    1, #Number of bands to create in the output
                    gdal.GDT_Float32)

#Copy across projection and transform details for the output
out.SetProjection(dataset.GetProjectionRef())
out.SetGeoTransform(in_dataset.GetGeoTransform()) 

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

#And write our processed data
out_band.WriteArray(out_data)
3
  • 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
    Commented Jun 14, 2012 at 18:53
  • That's right @JEquihua - most visualisation tools have the option of setting values transparent.
    – om_henners
    Commented Jun 16, 2012 at 1:18
  • 1
    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
    Commented Jun 18, 2012 at 4:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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