plot_corine_250.py

#!/usr/bin/env python
# conda create -n georasters 'georasters=0.5.15' ipython
import numpy as np
import georasters as gr
# from osgeo import gdalnumeric, gdal
from osgeo import gdal
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap


resize = False
extract = True
extract_palette = True

filename = '/home/megies/corine/2019-11-14/GLC2000_EU_250m.tif'
# ETRS89
extract_x = 4500000
extract_y = 2800000
extract_r = 50000


data = gr.from_file(filename)

# this definitly might need to be changed depending on data format of geotiff..
assert data.raster.dtype == np.uint8
vmin = np.iinfo(np.uint8).min
vmax = np.iinfo(np.uint8).max
# compare `np.iinfo(np.uint8)` and see `vmin`/`vmax` in ax.imshow below,
# otherwise the custom paletted colormap won't match the data


# set correct extent (might be wrong by half the pixle width at edges..)
xmin, dx, _, ymax, _, dy = data.geot
assert dy < 0  # dy usually is negative, if not needs np.flipud() on data
dy = -dy
ny, nx = data.raster.shape[-2:]
xmax = xmin + nx * dx
ymin = ymax - ny * dy
extent = (xmin, xmax, ymin, ymax)

# since .extract doesn't update geot, i guess resize won't do it either, so
# extent for plot needs to be computed before resize
# optionally reduce resolution, doesn't work to well when plotting the data
# paletted though it seems
if resize:
    old_shape = data.raster.shape
    new_shape = [int(i / 4.0) for i in old_shape]
    data = data.resize(new_shape)

# or extract a subraster
if extract:
    print(data.geot)
    data = data.extract(extract_x, extract_y, extract_r)
    # sadly geot doesn't get updated on the extract call
    print(data.geot)
    extent = [extract_x - extract_r, extract_x + extract_r,
              extract_y - extract_r, extract_y + extract_r]

# georasters does not take into account palette stored in geotiff..
data.plot()
#data.plot(extent=extent)

if extract_palette:
    # extract paletted colormap from geotiff with gdal
    # gdnum = gdalnumeric.LoadFile(filename)
    fn = gdal.Open(filename)
    band = fn.GetRasterBand(1)
    ct = band.GetRasterColorTable()
    num_colors = ct.GetCount()
    color_list = [ct.GetColorEntry(i) for i in range(num_colors)]
    # normalize RGBA values to 0-1
    # this might differ if palette is not stored as uint8 (0-255)..
    color_list = np.array(color_list) / 255.0
    cmap = LinearSegmentedColormap.from_list('geotiff', color_list, N=num_colors)
else:
    vmax = None
    vmin = None
    cmap = 'viridis'

fig, ax = plt.subplots()
ax.imshow(data.raster, cmap=cmap, vmin=vmin, vmax=vmax, extent=extent)
ax.grid()
plt.show()