#! /usr/bin/python3
# Last edited on 2019-05-09 01:33:51 by stolfilocal

# Read a GeoTIFF file, write as a 16-bit PGM file
#
# Required libraries:
#      numpy
#      rasterio
#      affine
#      matplotlib
#

import numpy as np
import rasterio
from affine import Affine
import matplotlib.pyplot as plt
import matplotlib.image as mping
import matplotlib.colors as colors

# directory path and file name
path = ''
file_in = 'GEBCO_2019_-7.0_65.0_38.0_25.0_Geotiff.tif'

# read in file with geotiff geographic information
src = rasterio.open(path + file_in)

# print out metadata information
for k in src.meta:
  print(k,src.meta[k])

# Retrieve the affine transformation
if isinstance(src.transform, Affine):
     transform = src.transform
else:
     transform = src.affine

N = src.width
M = src.height
dx = transform.a
dy = transform.e
minx = transform.c
maxy = transform.f

# Read the image data, flip upside down if necessary
data_in = src.read(1)
if dy < 0:
  dy = -dy
  data_in = np.flip(data_in, 0)

print('Data minimum, maximum = ', np.amin(data_in), np.amax(data_in))

print("P2\n")
print(N, " ", M, "\n");
print(2**16-1, "\n")

extent = [xdata[0], xdata[-1], ydata[0], ydata[-1]]
plt.figure(figsize=(8,8))
fig = plt.imshow(data_scale, extent=extent, origin='lower', cmap=cmap)
plt.title(file_in)
# Hide the axes and remove the space around them
plt.axis('off')
fig.axes.get_xaxis().set_visible(False)
fig.axes.get_yaxis().set_visible(False)
tickval = [1, 10, 100, 1000, 3000]
t = np.log(tickval)
cb = plt.colorbar(fig, ticks=255*(t - t[0])/(t[-1] - t[0]), shrink=0.5)
cb.set_label('Velocity Magnitude (m/year)')
cb.ax.set_yticklabels(tickval)
plt.savefig("GEBCO-py.png", dpi=300, bbox_inches='tight', pad_inches=0.5)
plt.show()