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
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]:
The dataset used here includes the following categories:
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()
In [41]:
bend = data[92000:98500,31500:35500]
plt.figure(figsize=(12,10))
bounds = np.linspace(0.5,7.5,8)
plt.imshow(bend, cmap=cmap,interpolation=None)
plt.plot([250,250+2*333],[6300,6300],'k',linewidth=3)
plt.text(250+140,6250,'20 km',color='k',fontsize=14)
plt.axis('off')
plt.colorbar(boundaries=bounds,shrink=0.3,aspect=20,orientation='horizontal',pad=0.01,ticks=[1,2,3,4,5,6,7])
plt.xlim(0,np.shape(bend)[1])
plt.ylim(np.shape(bend)[0],0)
plt.tight_layout()
In [55]:
csik = data[94100:96200,32500:35400]
plt.figure(figsize=(12,10))
bounds = np.linspace(0.5,7.5,8)
plt.imshow(csik, cmap=cmap,interpolation=None)
plt.plot([250,250+333],[2000,2000],'k',linewidth=3)
plt.text(250+100,1980,'10 km',color='k',fontsize=14)
plt.axis('off')
plt.colorbar(boundaries=bounds,shrink=0.3,aspect=20,orientation='horizontal',pad=0.01,ticks=[1,2,3,4,5,6,7])
plt.xlim(0,np.shape(csik)[1])
plt.ylim(np.shape(csik)[0],0)
plt.tight_layout()
In [56]:
csik_zoom = data[94850:95300,32700:33500]
plt.figure(figsize=(12,10))
bounds = np.linspace(0.5,7.5,8)
plt.imshow(csik_zoom, cmap=cmap,interpolation=None)
plt.plot([180,180+333/2.0],[430,430],'k',linewidth=3)
plt.text(180+60,425,'5 km',color='k',fontsize=14)
plt.axis('off')
plt.colorbar(boundaries=bounds,shrink=0.3,aspect=20,orientation='horizontal',pad=0.01,ticks=[1,2,3,4,5,6,7])
plt.xlim(0,np.shape(csik_zoom)[1])
plt.ylim(np.shape(csik_zoom)[0],0)
plt.tight_layout()