Introduction

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: http://glad.geog.umd.edu/europe/; it was put together by scientists at the University of Maryland -- for more detail, see this paper: http://dx.doi.org/10.1016/j.rse.2014.11.027. The notebook forms the basis of this blog post: http://hinderedsettling.com/2015/07/18/forest-cover-change-in-the-carpathians-since-1985/. 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))
    else:
        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]],
                         bit_rgb[colors[i][1]],
                         bit_rgb[colors[i][2]])
    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)
np.shape(data)
Out[4]:
(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)
plt.figure(figsize=(18,14))
bounds = np.linspace(0.5,7.5,8)
plt.imshow(carps, cmap=cmap,interpolation=None)
plt.colorbar(boundaries=bounds,shrink=0.3,aspect=20,orientation='horizontal',pad=0.01,ticks=[1,2,3,4,5,6,7])
plt.plot([250,250+666.7],[4300,4300],'k',linewidth=3)
plt.text(250+160,4250,'100 km',color='k',fontsize=14)
plt.axis('off')
plt.xlim(0,np.shape(carps)[1])
plt.ylim(np.shape(carps)[0],0)
plt.tight_layout()