import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from osgeo import gdal
import numpy as np

fig = plt.figure(figsize=(15,15))  # a new figure window
ax = fig.add_subplot(1, 1, 1)  # specify (nrows, ncols, axnum)
ax.set_title('Northern Norway', fontsize=14)

# Read the data and metadata
datafile = gdal.Open(r'2c_lavhei_68N_Jenny_wgs84.tif')
bnd1 = datafile.GetRasterBand(1).ReadAsArray()
nx = datafile.RasterXSize # Raster xsize
ny = datafile.RasterYSize # Raster ysize

gt = datafile.GetGeoTransform()
proj = datafile.GetProjection()

print("Geotransform",gt)
print("proj=", proj)
xres = gt[1]
yres = gt[5]

# get the edge coordinates and add half the resolution
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmax = gt[0] + (xres * nx) - xres * 0.5
ymin = gt[3] + (yres * ny) + yres * 0.5
ymax = gt[3] - yres * 0.5
print("xmin=", xmin,"xmax=", xmax,"ymin=",ymin, "ymax=", ymax)


# create a grid of lat/lon coordinates in the original projection
(lon_source,lat_source) = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
print(xmin,xmax+xres,xres, ymax+yres,ymin,yres)

i,j = np.unravel_index(bnd1.argmax(), bnd1.shape)
print(lon_source.shape, lat_source.shape)
print(lon_source[j,i], lat_source[j,i])

"""
map = Basemap(
    llcrnrlat=65.0,
    llcrnrlon=10.0,
    urcrnrlat=72.0,
    urcrnrlon=35.0,
    area_thresh=100.0,
    resolution='i',
    projection='lcc',
    lat_1=60., lon_0=15,
    rsphere=(6378137.00,6356752.3142))

# set missing values to 0 (for plotting...)
bnd1[bnd1 < 0] = 0
# project in the original Basemap and plot with pcolormesh
x,y = map(lon_source, lat_source)
map.pcolormesh(x,y,bnd1.T, cmap='spectral')
map.drawcountries(color='blue', linewidth=1.5, ax=ax)
map.drawcoastlines(linewidth=1.5, color='red', ax=ax)
plt.show()

"""