This is a notebook that I used to create maps from a Landsat-based dataset that documents forest cover change in Eastern Europe between 1985 and 2012. The dataset can be downloaded here:; it was put together by scientists at the University of Maryland -- for more detail, see this paper: The notebook forms the basis of this blog post: You need to have the GDAL package installed to run the analysis below.

In [57]:
import gdal
from gdalconst import *
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

%matplotlib inline

# function for making custom colormap (from Stackexchange)
def make_cmap(colors, position=None, bit=False):
    make_cmap takes a list of tuples which contain RGB values. The RGB
    values may either be in 8-bit [0 to 255] (in which bit must be set to
    True when called) or arithmetic [0 to 1] (default). make_cmap returns
    a cmap with equally spaced colors.
    Arrange your tuples so that the first color is the lowest value for the
    colorbar and the last is the highest.
    position contains values from 0 to 1 to dictate the location of each color.
    import matplotlib as mpl
    import numpy as np
    bit_rgb = np.linspace(0,1,256)
    if position == None:
        position = np.linspace(0,1,len(colors))
        if len(position) != len(colors):
            sys.exit("position length must be the same as colors")
        elif position[0] != 0 or position[-1] != 1:
            sys.exit("position must start with 0 and end with 1")
    if bit:
        for i in range(len(colors)):
            colors[i] = (bit_rgb[colors[i][0]],
    cdict = {'red':[], 'green':[], 'blue':[]}
    for pos, color in zip(position, colors):
        cdict['red'].append((pos, color[0], color[0]))
        cdict['green'].append((pos, color[1], color[1]))
        cdict['blue'].append((pos, color[2], color[2]))

    cmap = mpl.colors.LinearSegmentedColormap('my_colormap',cdict,256)
    return cmap

Load data

In [4]:
dataset = gdal.Open('/Users/zoltan/Downloads/dynamics_type/dynamics_type.img', GA_ReadOnly )
band = dataset.GetRasterBand(1)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
driver = dataset.GetDriver().LongName
data = band.ReadAsArray(0, 0, cols, rows)
(114000, 118000)

Create and display maps

The dataset used here includes the following categories:

  1. stable non-forest;
  2. stable forest;
  3. forest gain over non-forest in 1985;
  4. forest loss;
  5. forest loss followed by forest gain;
  6. repeated forest loss separated by forest gain;
  7. forest loss on areas which gain forest cover after non-forest state in 1985.

Whole Carpathian area

In [40]:
carps = data[79000:101500:5,13500:37000:5]
cmap = tuple([[0.85,0.85,0.85],[0.0,0.7,0.0],[0,0,1],[1,0,0],[0.3,0,0.8],[1,0,0.5],[0.6,0,0.2]])
cmap = make_cmap(cmap)
bounds = np.linspace(0.5,7.5,8)
plt.imshow(carps, cmap=cmap,interpolation=None)
plt.text(250+160,4250,'100 km',color='k',fontsize=14)