#!/usr/bin/env python
# coding: utf-8

# ### [WUR GRS-33806 Geoscripting](https://geoscripting-wur.github.io/)<img src='http://www.wageningenur.nl/upload/f9a0b5f2-15c5-4b84-9e14-20ef02f5e265_wur-logo.png' align="right">
# # Week 3 (Python): Raster handling with Python

# Good morning! Today we will handle some rasters in Python.

# ## Python modules for raster data handling

# The modules used below are:
# 
# * ***`gdal`***: bindings to the GDAL library which is part of the osgeo library
# * `numpy` for array calculations.
# 
# These two libraries are the base of all raster processing in Python!
# 
# There are additional libraries that offer wrapper functions for `gdal` and provide additional raster procssesing functionality:
# * mapbox [`rasterio`](https://github.com/mapbox/rasterio), see a [NDVI tutorial](http://www.loicdutrieux.com/pyLandsat/NDVI_calc.html)
# * [`rios`](http://rioshome.org/), a set of Python modules which makes it easy to write raster processing code in Python
# 
# Many important processing methods are available through additional libraries:
# * Python interface for Orfeo toolbox (`otb`), e.g. segmentation
# * `rsgislib` e.g. segmentation
# * `scikit-image` for image interpretation
# * `scikit-learn` for machine learning

# ## Data

# In the below example we will use an Aster image which you can download from the following [dropbox](https://www.dropbox.com/s/rsc4lzkd3t2adq5/ospy_data5.zip?dl=0). The data is also available [here](http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_data5.zip). 

# ## Reading, deriving NDVI, and writing raster data

# Let's get started and we will go through the following steps:
# 
# * Open the Aster image (ERDAS *.img format)
# * Read in the image data as an array
# * Derive the NDVI using array calculations
# * Write the resulting file as a *.tif (i.e. GeoTIFF file)
# * Close all files

# In[1]:


# import modules
from osgeo import gdal
from osgeo.gdalconst import GA_ReadOnly, GDT_Float32
import numpy as np

# open file and print info about the file
# the "./" refers to the parent directory of my working directory

filename = './data/ospy_data5/aster.img'
dataSource = gdal.Open(filename, GA_ReadOnly)

print("\nInformation about " + filename)
print("Driver: ", dataSource.GetDriver().ShortName,"/", \
      dataSource.GetDriver().LongName)
print("Size is ",dataSource.RasterXSize,"x",dataSource.RasterYSize, \
      'x',dataSource.RasterCount)

print('\nProjection is: ', dataSource.GetProjection())

print("\nInformation about the location of the image and the pixel size:")
geotransform = dataSource.GetGeoTransform()
if not geotransform is None:
    print('Origin = (',geotransform[0], ',',geotransform[3],')')
    print('Pixel Size = (',geotransform[1], ',',geotransform[5],')')


# In[2]:


# Read data into an array
band2Arr = dataSource.GetRasterBand(2).ReadAsArray(0,0,dataSource.RasterXSize, dataSource.RasterYSize)
band3Arr = dataSource.GetRasterBand(3).ReadAsArray(0,0,dataSource.RasterXSize, dataSource.RasterYSize)
print(type(band2Arr))
                                                   

# set the data type
band2Arr=band2Arr.astype(np.float32)
band3Arr=band3Arr.astype(np.float32)

# Derive the NDVI
mask = np.greater(band3Arr+band2Arr,0)

# set np.errstate to avoid warning of invalid values (i.e. NaN values) in the divide 
with np.errstate(invalid='ignore'):
    ndvi = np.choose(mask,(-99,(band3Arr-band2Arr)/(band3Arr+band2Arr)))
print("NDVI min and max values", ndvi.min(), ndvi.max())
# Check the real minimum value
print(ndvi[ndvi>-99].min())


# In[3]:


# Write the result to disk
driver = gdal.GetDriverByName('GTiff')
outDataSet=driver.Create('./data/ndvi.tif', dataSource.RasterXSize, dataSource.RasterYSize, 1, GDT_Float32)
outBand = outDataSet.GetRasterBand(1)
outBand.WriteArray(ndvi,0,0)
outBand.SetNoDataValue(-99)

# set the projection and extent information of the dataset
outDataSet.SetProjection(dataSource.GetProjection())
outDataSet.SetGeoTransform(dataSource.GetGeoTransform())

# Finally let's save it... or like in the OGR example flush it
outBand.FlushCache()
outDataSet.FlushCache()


# Let's check the resulting file via *Bash* (all commands in the current document with a **!** in front are run via the system shell, which is *Bash*). Below we call the `gdalinfo` command of the **GDAL** C++ library directly. If you want to learn more about `gdalinfo` go to: http://www.gdal.org/gdalinfo.html

# In[4]:


get_ipython().system('gdalinfo ./data/ndvi.tif')


# Important here is that the `gdalinfo` contains information about the projection  (via `SetProjection`) and also about the corner coordinates (via `SetGeoTransform`).

# ## Reproject the NDVI image to Lat/long

# Reprojecting a raster image can be done by
# 
# * calling the `gdalwarp` function of the **GDAL C++** library directly
# * using `Python` and the `gdal` module
# 
# ### Using `gdalwarp`
# 
# The easiest way to reproject a raster file is to use GDAL's [gdalwarp](http://www.gdal.org/gdalwarp.html) tool. As an example, we will reproject the above NDVI image derived earlier into latitude/longitude (WGS84).
# 
# * `-t_srs "EPSG:4326"` is the CRS to reproject **to**, i.e. lat/long, known by its [EPSG code](http://spatialreference.org/ref/epsg/4326/). The CRS to reproject **from** is already specified in the source data set, see above.
# * `data/ndvi.tif` is the **input** dataset.
# * `data/ndvi_ll.ti`f is the **output** dataset (the output format is here automatically set by the extension i.e. GeoTIFF).
# 
# > **Question 1**: What is the CRS of the image before it is reprojected?

# In[7]:


# via the Shell
get_ipython().system('gdalwarp -t_srs "EPSG:4326" ./data/ndvi.tif ./data/ndvi_ll.tif')


# In[8]:


# Let's check what the result is
get_ipython().system('gdalinfo ./data/ndvi_ll.tif')


# To reproject a raster from within Python is not as straight-forward... Have a look at the following links, and use the internet to find a way!
# 
# - The function defined by Prof. P. Lewis (see here [for more information](http://nbviewer.ipython.org/github/profLewis/geogg122/blob/master/Chapter4_GDAL/GDAL_Python_bindings.ipynb)). 
# - Core GDALDataset Class reference info (see [GeoTransFrom Info](http://www.gdal.org/classGDALDataset.html#a0fe0f81d65d84557b5d71ddc024faa02). For a relevant question and example on GIS Stack Exchange ([see](https://gis.stackexchange.com/questions/24055/raster-data-array-output-flipped-on-x-axis-using-python-gdal)).
# - [GIS Stack Exchange](https://gis.stackexchange.com/questions/6669/converting-projected-geotiff-to-wgs84-with-gdal-and-python).
# - Or use other raster libraries, e.g. [rasterio](https://github.com/mapbox/rasterio/blob/master/examples/reproject.py).

# ### Visualise the result
# 
# Let's see the output of the reprojected NDVI image.

# In[9]:


# Notebook magic to select the plotting method
# Change to inline to plot within this notebook
#%matplotlib inline 
get_ipython().run_line_magic('matplotlib', 'inline')
from osgeo import gdal
import matplotlib.pyplot as plt


# In[11]:


# Open image
dsll = gdal.Open("./data/ndvi_ll.tif")

# Read raster data
ndvi = dsll.ReadAsArray(0, 0, dsll.RasterXSize, dsll.RasterYSize)

# Now plot the raster data using gist_earth palette
plt.imshow(ndvi, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth)
plt.show()

dsll = None


# ### Check in `R`
# 
# We can quickly check what we have done with R (go now to your Rstudio, or R terminal) and try the following code:
# 
# ```
# library(raster)
# a <- brick('./data/ospy_data5/aster.img')
# plotRGB(a, 3, 2, 1, stretch='hist')
# ```
# 
# ```
# b <- raster('./data/ndvi.tif')
# extent(b)
# projection(b)
# hist(b, 1, maxpixels=500, plot=TRUE)
# ```

# ## Exercise lesson 12: NDWI and reprojection in python
# 
# For this scripting challenge I have downloaded a Landsat image with all bands processed to surface reflectance (Level 1T). You can download it from [here](https://www.dropbox.com/s/zb7nrla6fqi1mq4/LC81980242014260-SC20150123044700.tar.gz?dl=0). Unzip the file and you will see that it contains all the invidual bands. 
# 
# Now, write a well-structured and documented script where you define ***a function to derive the Normalised difference water index (NDWI)*** and derive it from the landsat image and ***reproject the image to Lat/Long WGS84*** (**hint**: use GDAL from *Bash* in your notebook).
# 
# `NDWI = band 3 - band 5 / band 3 + band 5` ([more info about](http://nsidc.org/data/docs/daac/nsidc0332_smex03_srs_ndvi_ndwi_ok.gd.html) `NDWI`)
# 
# A clean and well structured script is critical here.
# 
# *** Create an Ipython Notebook file (.ipynb) with Jupyter Notebook and put it on GitLab!***
# 
# 

# ## Additional Resources
# - Chris Holden from Boston University on [handling raster data with GDAL](http://www.gis.usu.edu/~chrisg/python/2009/)
# - [Blog about Python for GeoSpatial data analysis](http://www.digital-geography.com/python-for-geospatial-data-analysis-part-ii/?subscribe=invalid_email#477)
# - https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html
# - Great basic tutorial: http://www.gdal.org/gdal_tutorial.html
# - [GDAL Python info and bindings](http://gdal.org/python/osgeo.gdal_array-module.html#BandReadAsArray)
# - yet another example https://borealperspectives.wordpress.com/2014/01/16/data-type-mapping-when-using-pythongdal-to-write-numpy-arrays-to-geotiff/
# - Typical problems with GDAL and python: [PythonGotchas](http://trac.osgeo.org/gdal/wiki/PythonGotchas)
# * http://www.qgistutorials.com/nl/docs/getting_started_with_pyqgis.html
# * [QGIS and Python tutorial](http://www.qgisworkshop.org/html/workshop/python_in_qgis_tutorial2.html)
# 

# In[ ]: