Plotting African ecosystems

data_dir = 'data/geo/af_labeled_ecosys'
!ls -lh {data_dir}

# conda install -c conda-forge gdal
from osgeo import osr, gdal
print('gdal', gdal.__version__)

# open the dataset
import os
dataset = gdal.Open(os.path.join(data_dir, 'Africa_IVC_20130316_final_MG.tif'))

proj_wkt = dataset.GetProjection()
print(proj_wkt)

proj = osr.SpatialReference()
proj.ImportFromWkt(proj_wkt)
print(proj)

geo_transform = dataset.GetGeoTransform()
geo_transform

origin_x = geo_transform[0]
origin_y = geo_transform[3]
pixel_width = geo_transform[1]
pixel_height = geo_transform[5]

# how big?
n_cols = dataset.RasterXSize
n_rows = dataset.RasterYSize
n_cols, n_rows

# how many bands?
dataset.RasterCount

band = dataset.GetRasterBand(1)

# what data type?
gdal.GetDataTypeName(band.DataType)

import humanize
humanize.naturalsize(n_cols * n_rows * 4)

band.GetNoDataValue()

band.GetMinimum()

band.GetMaximum()

band.GetColorTable() is None  # sadly

# apt-get install dbview
!dbview {data_dir}/Africa_IVC_20130316_final_MG.tif.vat.dbf --browse --trim --description > \
    {data_dir}/Africa_IVC_20130316_final_MG.tif.vat.csv

!head -n23 {data_dir}/Africa_IVC_20130316_final_MG.tif.vat.csv

# extract the field names
import petl as etl
vat_fn = os.path.join(data_dir, 'Africa_IVC_20130316_final_MG.tif.vat.csv')
tbl_descr = etl.fromtsv(vat_fn).head(19).convertall('strip')
tbl_descr.displayall()

# extract the data
hdr_colors = tbl_descr.values('Field Name').list()
tbl_colors = (etl
    .fromcsv(vat_fn, delimiter=':')
    .skip(20)  # skip the field descriptions
    .pushheader(hdr_colors)
    .cat()  # remove empty cells beyond columns
    .convert('Value', int)
    .convert(['Count', 'Red', 'Green', 'Blue', 'Opacity'], float)
)
tbl_colors

import numpy as np
max_class = tbl_colors.values('Value').max()
max_class

colors = np.zeros((max_class+1, 3), dtype=float)
for i, r, g, b in tbl_colors.cut('Value', 'Red', 'Green', 'Blue').data():
    colors[i] = r, g, b

colors

# check all numbers processed ok
np.count_nonzero(np.isnan(colors))

# fix the zero colour to be white
colors[0] = 1, 1, 1

import matplotlib as mpl
color_map = mpl.colors.ListedColormap(colors)

data = dataset.ReadAsArray(buf_xsize=n_cols//100, buf_ysize=n_rows//100, 
                           resample_alg=gdal.GRIORA_Mode)
data

data.shape

import matplotlib.pyplot as plt
%matplotlib inline

# set all missing data as 0 class
data[data < 0] = 0

plt.imshow(data, cmap=color_map);

# conda install -c conda-forge cartopy
import cartopy; print('cartopy', cartopy.__version__)
import cartopy.crs as ccrs
import cartopy.feature as cfeature

extent_lonlat = (
    origin_x, 
    origin_x + (pixel_width * dataset.RasterXSize),
    origin_y + (pixel_height * dataset.RasterYSize),
    origin_y
)

print(proj)

crs_lonlat = ccrs.PlateCarree()

subplot_kw = dict(projection=crs_lonlat)
fig, ax = plt.subplots(figsize=(9, 9), subplot_kw=subplot_kw)
ax.set_extent(extent_lonlat, crs=crs_lonlat)
ax.imshow(data, cmap=color_map, extent=extent_lonlat, origin='upper')
ax.coastlines(resolution='50m', linewidth=.5)
gl = ax.gridlines(crs=crs_lonlat, 
                  xlocs=np.arange(-180, 180, 15), 
                  ylocs=np.arange(-180, 180, 15), 
                  draw_labels=True)
gl.xlabels_top = None
ax.add_feature(cfeature.BORDERS, linewidth=.5)
ax.set_title('Plate carrée projection', va='bottom');

crs_laea = ccrs.LambertAzimuthalEqualArea()
xmin, ymin = crs_laea.transform_point(extent_lonlat[0], extent_lonlat[2], crs_lonlat)
xmax, ymax = crs_laea.transform_point(extent_lonlat[1], extent_lonlat[3], crs_lonlat)
extent_laea = (xmin, xmax, ymin, ymax)

subplot_kw = dict(projection=crs_laea)
fig, ax = plt.subplots(figsize=(9, 9), subplot_kw=subplot_kw)
ax.set_extent(extent_laea, crs=crs_laea)
ax.imshow(data, cmap=color_map, extent=extent_lonlat, origin='upper', transform=crs_lonlat)
ax.coastlines(resolution='50m', linewidth=.5)
gl = ax.gridlines(crs=crs_lonlat, 
                  xlocs=np.arange(-180, 180, 15), 
                  ylocs=np.arange(-180, 180, 15))
ax.add_feature(cfeature.BORDERS, linewidth=.5)
ax.set_title('Lambert azimuthal equal area projection');

# conda install -c conda-forge basemap
# conda install -c conda-forge basemap-data-hires
from mpl_toolkits.basemap import Basemap
import mpl_toolkits.basemap; print('basemap', mpl_toolkits.basemap.__version__)

plt.figure(figsize=(9, 9))

# plate carrée projection is known as 'cyl' in basemap
map = Basemap(projection='cyl', resolution='l',
              llcrnrlon=extent_lonlat[0], llcrnrlat=extent_lonlat[2],
              urcrnrlon=extent_lonlat[1], urcrnrlat=extent_lonlat[3])

# plot the raster data
x = np.linspace(extent_lonlat[0], extent_lonlat[1], data.shape[1])
y = np.linspace(extent_lonlat[3], extent_lonlat[2], data.shape[0])
xx, yy = np.meshgrid(x, y)
map.pcolormesh(xx, yy, data, cmap=color_map)

# plot coastlines etc.
map.drawcoastlines(linewidth=.5)
map.drawcountries(linewidth=.5)
parallels = np.arange(-45, 45, 15)
map.drawparallels(parallels, labels=[1, 0, 0, 0])
meridians = np.arange(0, 360, 15)
map.drawmeridians(meridians, labels=[0, 0, 0, 1])
plt.title('Plate carrée projection', va='bottom');

plt.figure(figsize=(9, 9))

# use Lambert azimuthal equal area projection
m = Basemap(projection='laea', resolution='l', 
            lon_0=(extent_lonlat[0] + extent_lonlat[1])/2,
            lat_0=(extent_lonlat[2] + extent_lonlat[3])/2,
            llcrnrlon=extent_lonlat[0], llcrnrlat=extent_lonlat[2],
            urcrnrlon=extent_lonlat[1], urcrnrlat=extent_lonlat[3])

# plot the raster data, transforming coordinates using the map's projection
lons = np.linspace(extent_lonlat[0], extent_lonlat[1], data.shape[1])
lats = np.linspace(extent_lonlat[3], extent_lonlat[2], data.shape[0])
lons, lats = np.meshgrid(lons, lats)
xx, yy = m(lons, lats)
m.pcolormesh(xx, yy, data, cmap=color_map)

# plot coastlines etc.
m.drawcoastlines(linewidth=.5)
m.drawcountries(linewidth=.5)
parallels = np.arange(-45, 45, 15)
m.drawparallels(parallels, labels=[1, 0, 0, 0])
meridians = np.arange(0, 360, 15)
m.drawmeridians(meridians, labels=[0, 0, 0, 1])
plt.title('Lambert azimuthal equal area projection', va='bottom');