Editing .tiff rasters in Python

Question

I am trying to edit .tif rasters in python, but I have some issues regarding the tiff to array conversion. I read there are multiple options for opening rasters (gdal, PIL, matplotlib), but which of those is the best for allowing me to convert to array, edit (apply local filters), then convert back to .tif?

For example, I get a strange result when opening arrays from a .tif raster using PIL library.

The code looks like this:

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import matplotlib.cm as cm
import numpy as np
import pywt as pw
import tkFileDialog as tk
from PIL import Image

def fct(fisier):
    img = Image.open(fisier)
    arr=np.array(img, dtype=np.float)
    plt.imshow(arr, cmap=cm.Greys_r)
    plt.show()
    print(arr)
    return arr, img
file = tk.askopenfile(initialdir='C:/temp', mode='rb')
fct(file)

But the array that is displayd in all white, and has this form:

[[ -2.14748365e+09  -2.14748365e+09  -2.14748365e+09 ...,   0.00000000e+00
    0.00000000e+00   6.19762981e+08]

What do I need to do in order to open the tiff raster as an array, so I can edit it further?

Answer 1

Image libraries such as PIL or opencv are great for image processing but ignore the geolocation side of geotiffs. We want to use gdal because it recognizes most of the projection types and it is able to save them within a geotiff file that will be compatible with other GIS software (for example, both ArcGIS and QGIS use python gdal).

Several modules that handle geographic tiff files (i.e. a raster with metadata containing geolocation) are definitely more user friendly and intuitive than gdal. However, most (if not all) of them USE gdal.

When it comes to modify pixel values after some calculation (say, divide every pixel of a geotiff raster by 4), my workflow can be summarized as below:

1. Open the tiff with gdal
2. Read and store size, geotransform and projection
3. Load the raster as numpy array
4. Execute calculation with numpy
5. Create an empty tiff container with same size, geotransform and projection of original
6. Save the numpy array in the tiff file

Note: geotransform is an object that stores pixel size and coordinates of bottom left corner pixel according to the reference system stored in the projection object.

import gdal

#1.    
tiff_file = gdal.Open(raster_name)
#2.
geotransform = tiff_file.GetGeoTransform()
projection = tiff_file.GetProjection()
band = tiff_file.GetRasterBand(1)    
xsize = band.XSize
ysize = band.YSize
#3.
array = band.ReadAsArray()
tiff_file = None #close it
band = None #close it
#4.
array = array/4
#5.
driver = gdal.GetDriverByName('GTiff')
new_tiff = driver.Create(output_raster_name,xsize,ysize,1,gdal.GDT_Int16)
new_tiff.SetGeoTransform(geotransform)
new_tiff.SetProjection(projection)
new_tiff.GetRasterBand(1).WriteArray(array)
new_tiff.FlushCache() #Saves to disk 
new_tiff = None #closes the file

Even though the gdal library might not have the greatest documentation (in terms of readability for python users), the methods have pretty straightforward naming. This can be of great help once you know what you want to do.

Answer 2

With GDAL:

from osgeo import gdal
file = gdal.Open(’raster1.tif’)
band1 = file.GetRasterBand(1)
a = band1.ReadAsArray()

In tiff you have 3 bands (RGB). So after you open raster file, you have to get band (1,2 or 3 are Red, Green or Blue) and then read band as array. Of course there could be only one band, but even then you have to get it before readAsArray.

Answer 3

GDAL is the most widely used library to access geospatial data from python. Unfortunately it can be a bit unintuitive from time to time. A recent development trying to make the use of GDAL more pythonic is rasterio

import rasterio

with rasterio.open('path/to/your/geo.tiff') as ds:
        ds = src.read()