plot_geotiff.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

filename = 'PRINT_DTK25_Rosenheim.tif'
resize = False
extract = True

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
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)
    extract_x = 4507000
    extract_y = 5302000
    extract_r = 1000
    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(extent=extent)

# 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.
cmap = LinearSegmentedColormap.from_list('geotiff', color_list, N=num_colors)

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