%matplotlib inline

"""
Explanation: Object based training and classification using machine learning

A simple land cover mapping example using aerial photography
The following workflow demonstrates using image segmentation and machine learning to produce landcover classes, the priniciples of which broadly apply to much remote sensing work. Image segmentation  then classification of the resulting attributes is often used to reduce noisy/detailed imagery to lancover classes greater in spatial coverage than the detail within the image.   
The data:


Colour infra-red imagery collected over ther Arun Valley, Sussex (G, R, NiR), courtesy of the now defunct LandMap website


shapefiles of the training samples


Disclaimer! I do not claim this is good model or set of classes!! 
This is merely to demonstrate the funtionality on a segmentation derived from fine spatial resolution data
Obviously, the reality of this type of approach is using large training sets and lengthier model parameter searches are required, as datasets are bigger than this small demonstration. 
The Data
The data is available here
https://drive.google.com/file/d/12GFv8473HdE1KoW0ESFXZgXCCGFysP0A/view?usp=sharing
The zip contains:

Arundel_WWT.tif (a Green/Red/NIR composite)
ArunLcover.qml (a QGIS colour map for the results)
Arundel_seg_empty.shp (the result of segmentation performed in this analysis complete with training labels)
End of explanation
"""


from geospatial_learn import learning, shape, raster
import shapefile

cd 'my/path/LCover'

"""
Explanation: Module imports

As mentioned above, the code above imports the modules required for this image classification workflow, using the Object-based image analysis method on polygonal data.
We use pyshp breifly to look at some attributes
End of explanation
"""


msRas = ('Arundel_WWT.tif')

segShape = ('Arundel_seg_empty.shp')


"""
Explanation: Input file paths
End of explanation
"""


# zonal stats
# please note that by using enumerate we assume the bandnames are ordered as the are in the image!
bandnames = ['g', 'r', 'nir']


# Please note we add 1 to the bnd index as python counts from zero
for bnd,name in enumerate(bandnames):
    shape.zonal_stats(segShape, msRas, bnd+1, name+'mn', stat = 'mean', write_stat = True)
    shape.zonal_stats(segShape, msRas, bnd+1, name+'mdn', stat = 'median', write_stat = True)
    shape.zonal_stats(segShape, msRas, bnd+1, name+'skw', stat = 'skew', write_stat = True)
    shape.zonal_stats(segShape, msRas, bnd+1, name+'krt', stat = 'kurtosis', write_stat = True)
    
   

# texture props
bandnames = ['g', 'r', 'nir']

for bnd,name in enumerate(bandnames):
    shape.texture_stats(segShape, msRas,  bnd+1,  gprop='entropy',  write_stat=True)



# shape props

props = ['MajorAxisLength', 'MinorAxisLength', 'Area', 'Eccentricity', 'Solidity',
         'Extent', 'Perimeter']
for prop in props:
    shape.shape_props(segShape, prop, label_field = 'DN')

"""
Explanation: Training & Imagery
The training and imagery should look something like this by setting importing the symbology from ArunLcover.qml.
Remeber to set the attribute to Train as this qml is set to RF, which is for the later results!
<img src="figures/ArunTrain.png" style="height:400px">
The classes used are:
<img src="figures/ArunClasses.png" style="height:100px">
The segmentations that was used
Included for reference only, I don't claim this algorithm/lib is the best solution, just an example
The example segmentation uses the OTB meanshift. 
The otb command for reference
bash
otbcli_Segmentation -in path/to/image -filter meanshift -filter.meanshift.spatialr 5 -filter.meanshift.ranger 10 
-filter.meanshift.minsize 50 -mode vector -mode.vector.out path/to/outShape                     
We prefix the command line with ! to execute in this notebook as this is external to python
Only run this if you wish to either experiment with parameters and/or wish to label the training segments yourself! Uncomment if this is the case.
Collection of statistics

Segment attribute data is written to the labeled shapefile using the shape.zonal_stats and shape.texture_stats function from geospatial_learn.
Rather than repeat the same line of code for each band, a simple for loop is used below to extract the band data.
End of explanation
"""


r = shapefile.Reader(segShape)

feat_fields = [f[0] for f in r.fields]
feat_fields = feat_fields[3:25]
feat_fields


len(feat_fields)

"""
Explanation: Model training
Now that we have collected our training statistics, we can calibrate our chosen model, which, in this case is a random forest.
Geospatial_learn wraps scikit-learn (hence the name) and xgboost, two excellent machine learning libraries (I intend to add tensor flow to this shortly also!).
Typically, we collect training from a shapefile and save it as a .gz file with the get_training_shp function, though it is possible to feed the training array straight into the create_model function.
The label and feature fields respectively. 
```python
label_field = 'Train'
feat_fields = ['gmn','gmdn','gskw','gkrt','rmn','rmdn','rskw','rkrt','nirmn','nirmdn',
                  'nirskw','nirkrt','entropy','entropy_1','entropy_2']
```
Path to the training array we intend to save. This is optional, it is possible to run
```python
training = 'savemytraining.gz'
```
```python
dfTrain, rejects = learning.get_training_shp(segShape, label_field, feat_fields,  outFile = training)
```
Motivational lazyness
Here is a quick way of getting a list of the fields to be used as features in the training and classification without typing them out
Here a list comprehension is used with the pyshp (shapefile) lib. This is not really that important, but if you are interested, it is a useful python concept....
If you are not familiar with list comprehensions, here is an attempt at explanation in the context of the shapefile. List comprehensions are useful, comapct and readable way of performing an operation on data.
The shapefile is read in as a python class object. A class is a fundamental part of object-based programming  (nothing to with the current task!!!) where the data structure has: 


properties - effectively attributes of the data


methods - built in operations on the data


Below r is the 'object' read in by the pyshp library
```python
r = shapefile.Reader(segShape)
The shapefile has a 'fields' property, but this returns a list in which each entry has the field name, datatype length, value. So each entry contains superfluous info:python
["AREA", "N", 18, 5],
```
Where as we are only interested in the first part of every entry, which is the field name. 
so this code is basically saying, output the first part (f[0]) of each entry (f) in the fields property (r.fields)
python
[f[0] for f in r.fields]
Hence we end up with a list of only the field titles. Which gives us the names of all the properties.
Finally we exclude the DN and Train fields by indexing.
python
feat_fields[3:25]
End of explanation
"""


training = 'Arundel_train_ob.gz'

label_field = 'Train'

dfTrain, rejects = learning.get_training_shp(segShape, label_field, feat_fields,  outFile = training)

"""
Explanation: Now to create some training samples
These can be used in both a pixel and object based manner
The variables returned are a dataframe (like in R) and a list of reject polygons which may have had invalid geometry.
In practice there shouldn't be many of these - it's just to make sure!
End of explanation
"""


dfTrain.head()

"""
Explanation: Here's a quick look at the returned dataframe to see the information. Dataframes are handy to visualise in jupyter and ipython. The create_model accepts the raw array data though so the .as_matrix() property of the dataframe is used when inputting to create_model
The first column is the label, the rest features
End of explanation
"""


params = {'n_estimators': [500], 'max_features': ['sqrt', 'log2'], 
          'min_samples_split':[5,10,20,50], 'min_samples_leaf': [5,10,20,50]}

outModel = 'Arundel_Rf_model2.gz'

"""
Explanation: Model creation with k-fold cross validated grid search
Now we can train the model with the above data.
We first define the parameters we wish to grid search over. The parameters below are just an example, It is of course possible for these to be more numerous at the cost of processing time. The time is a function of the number of possibilities per parameter.
```python
params = {'n_estimators': [500], 'max_features': ['sqrt', 'log2'], 
          'min_samples_split':[5,10,20,50], 'min_samples_leaf': [5,10,20,50]}

```        
When we execute the create_model function we get a summary of the no of model fits
'Fitting 5 folds for each of 18 candidates, totalling 90 fits'
I have fixed the n_estimators (trees) at 500 below but this could be varied also.
For a full list of params see:
http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html
End of explanation
"""


learning.create_model(dfTrain.as_matrix(), outModel, clf='rf', cv=5, params=params)

"""
Explanation: The create_model function is executed below. The progress of the gird search is printed out.
End of explanation
"""


learning.classify_object(outModel, segShape, feat_fields, field_name='RF')

"""
Explanation: At this point it is worth checking the shapefile to see everything is written correctly in qgis.....
Classification
Having collected training and created a model, we can finally classify our polygon attributes.
This is done using the classify object function where:
```python
learning.classify_object(outModel, inShape, feat_fields, field_name='RF')
```
We reuse the outModel, inShape and feat_fields variables from earlier which leaves only the keyword arg field_name, which is what we intend to call the field holding the classification values.
Keep this short if writing to ESRI shapefiles as there is a strict limit to the no of characters.
End of explanation
"""


learning.plot_feature_importances(outModel, feat_fields)

"""
Explanation: Model evaluation
Have a look at the results in QGIS by using the style file provided
<img src="figures/ArunRF.png" style="height:400px">
<img src="figures/ArunClasses.png" style="height:100px">
Having classified the segmentation - it looks fairly convincing, though some things appear strange, such as slivers of area in the fields classified as trees - could be explored with feature importance. 
The classifier may be splitting nodes on the basis of a geometric feature which is not actually indicative of trees!
Additionally part of the ponds on the top right are mistaken as impervious, though  superficially they are a similar colour.
End of explanation
"""


