import matplotlib.pyplot as plt
import numpy as np

import cartopy.crs as ccrs
from cartopy.io.img_tiles import GoogleTiles
from cartopy.io.srtm import srtm_composite

from osgeo import gdal
from osgeo import gdal_array

# Specify a region of interest, in this case, Sudelfeld Ski Resort (Germany)
lat = 47 + 40 / 60.0 + 30 / 3600.
lon = 12 + 3 / 60.0 + 2 / 3600.

plt.figure(figsize=(10, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.set_extent([12.0, 13.0, 47.0, 48.0])
gg_tiles = GoogleTiles()
ax.add_image(gg_tiles, 10)

plt.scatter(lon, lat, marker=(5, 1), color='red', s=200)
plt.title("Welcome to Sudelfeld")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False

plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())

elev, crs, extent = srtm_composite(12, 47, 1, 1)
plt.imshow(elev, extent=extent, transform=crs,
           cmap='gist_earth', origin='lower')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("SRTM Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False

plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())

elev, crs, extent = srtm_composite(12, 47, 1, 1)
elev = np.ma.masked_less_equal(elev,0,copy=False)
plt.imshow(elev, extent=extent, transform=crs,
           cmap='gist_earth', origin='lower')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("SRTM Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False

plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())

elev, crs, extent = srtm_composite(12, 47, 1, 1)
src_ds = gdal_array.OpenArray(elev)
srcband = src_ds.GetRasterBand(1)
dstband = srcband
maskband = srcband
smoothing_iterations = 0
options = []
max_distance = 200
result = gdal.FillNodata(dstband, maskband,
                         max_distance, smoothing_iterations, options,
                         callback=None)
elev = dstband.ReadAsArray()

plt.imshow(elev, extent=extent, transform=crs,
           cmap='gist_earth', origin='lower')
cb = plt.colorbar(orientation='vertical')
cb.set_label('Altitude')
plt.title("SRTM Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False

## Subplot 23: Shaded Relief Map
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
x, y = np.gradient(elev)

slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))

# -x here because of pixel orders in the SRTM tile
aspect = np.arctan2(-x, y)

altitude = np.pi / 6.
azimuth = np.pi /2.

shaded = np.sin(altitude) * np.sin(slope)\
    + np.cos(altitude) * np.cos(slope)\
    * np.cos((azimuth - np.pi/2.) - aspect)

plt.imshow(shaded, extent=extent, transform=crs,
           cmap='Greys', origin='lower')

plt.title("Shaded Relief Map")
gl = ax.gridlines(draw_labels=True,)
gl.xlabels_top = False
gl.ylabels_left = False