Skip to content

Instantly share code, notes, and snippets.

@megies
Created November 14, 2019 11:20
Show Gist options
  • Save megies/9b7f31f617a1a4c077e83d013d27d1f9 to your computer and use it in GitHub Desktop.
Save megies/9b7f31f617a1a4c077e83d013d27d1f9 to your computer and use it in GitHub Desktop.
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment