bodem = np.load("../data/nete_bodem.npy") # Should be applied to a real DEM data set...

def shaded_relief(elev, altitude=np.pi/4.,
                  azimuth=np.pi/2.):
    """
    Make a shaded relief version of a numpy array.
    Parameters
    ----------
    dem : 2D numpy array
    altitude : float
        Altitude of sun (in radians)
    azimuth : float
        Direction of sun (in radians)
    """
    x, y = np.gradient(elev)

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

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

plt.imshow(shaded_relief(bodem))

"""
Explanation: Gdal provide many functionalities, but...
For a number of examples: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html
(indeed, not the sweetest implementations)
But it does not align well with rasterio...
<div class="alert alert-danger">
    <b>REMEMBER</b>: 
    <ul>
        <li> Do NOT combine `import gdal` and `import rasterio` in the same script/notebook/...
        <li> calling GDAL with a subprocess from the command line is not an issue
    </ul>
</div>

Hence, we do not include the gdal example here...
Numpy (if you can...)
Numpy and scipy are great resources!
Furthermore, lots of useful functions are available for specific actions. For example:
End of explanation
"""


import rasterio

with rasterio.open("../scratch/NE1_50M_SR/NE1_50M_SR_Oost_Vlaanderen.tif") as oostvl:
    img = oostvl.read(1)

plt.imshow(img)

"""
Explanation: For a practical example, check the draft exploration  of the Niche application submodules trofie and pH. A more advanced application with Satellite Images Classification as topic is provided in this tutorial.
rasterio is the pythonic implementation specific for rasters
Reading raster files
End of explanation
"""


type(img)

"""
Explanation: Again, it easily interacts with numpy arrays:
End of explanation
"""


import rasterio
from rasterio.features import rasterize
from geopandas import read_file
import os
import matplotlib.pyplot as plt

data_dir = os.path.join('..', 'scratch', 'data_scipy2015')

raster_file = os.path.join(data_dir, 'manhattan.tif')
vector_file = os.path.join(data_dir, 'nybb_15b', 'nybb.shp')

df = read_file(vector_file)
plt.figure(1)


def show_img(img, bounds):
    left, bottom, right, top = bounds
    plt.imshow(img, cmap='gray', extent=(left, right, bottom, top))
    locs, labels = plt.xticks()
    plt.setp(labels, rotation=90)

with rasterio.open(raster_file) as src:
    # extract Manhattan, transform coordinates and geometry only
    poly = df.to_crs(src.crs).ix[3]['geometry']
    coords = [p.exterior.coords.xy for p in poly]
    ax = plt.subplot(1, 2, 1)
    img = src.read(1)
    show_img(src.read(1), src.bounds)
    for x, y in coords:
        plt.plot(x, y, color='red', linewidth=2)
    ax.set_xlim(src.bounds.left, src.bounds.right)
    ax.set_ylim(src.bounds.bottom, src.bounds.top)
    plt.subplot(1, 2, 2)
    mask = rasterize([poly], transform=src.transform, out_shape=src.shape)
    # If we had opened the file with mode 'r!', we could write the mask
    # src.write_mask(mask)
    img[mask==0] = 255
    show_img(img, src.bounds)

plt.savefig('rasterize_mask.png', dpi=300)
plt.figure(2)
show_img(mask, src.bounds)
plt.savefig('mask.png', dpi=300)

"""
Explanation: Example: rasterize vector feature
(Remark, you have to download the example data first and put it in the scratch folder to use this script)
End of explanation
"""


