raster - Creating a mask with transparent pixels in Python

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.

Answer

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)