<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Biomass Estimation - Putting the RAT into TPOT
<img src='http
Step1:  Prepare the training data
This section reads in the Biomass library csv file and links the observations to the closest segment.
NOTE
Step2:  Build the Machine Learning Model
This uses the amazing TPOT tool to find a good machine learning pipeline that predicts biomass. It automates the most tedious part of machine learning by intelligently exploring thousands of possible pipelines to find the best one for your data. There are a few setting below to tweak.
  - optTime is the number of minutes you are prepared to wait for the optimisation. The longer the better.
  - nCPUs is the number of CPUs you have available for the search. You'll also need a minimum of about 2GB of RAM per CPU
  - Scoring is the thing you'll be minimising
  - nSubsets picks a random subset from the training data. Smaller is faster but you might miss some variation in the data
  - The field from the biomass library you want to fit. Here we're selecting the Total Biomass column (number 10)
This will run a genetic algorithm until it times out and selects the best fitting model to use for biomass prediction. It saves the model so you can play with it later if you want.
Step4:  Predict Biomass using the RAT
Once we have the biomassModel we can use that to predict the Biomass for every entry in the RAT. 
This is cool, because for Australia we only need to make around 3,000,000 predictions which is far faster than predicting every one of the 18,407,102,543 pixels in the image. That's around a 6000x speedup.
Here,  RIOS ratapplier does the work writing a brand new attribute column into the kea file.
Step5:  Export the Biomass band to an individual GeoTIFF
You can run this section if you want to output just the newly calculated biomass column to a separate GeoTIFF. 
This uses one of the function from the rastergis module in RSGISLIB
After this, you might like to run gdal_translate to compress the output and gdaladdo to add overviews for easy storage and fast viewing in your favorite GIS platform. Like this
<ASSISTANT_TASK:>
Python Code:
# Imports for this Python3 notebook
import numpy
import matplotlib.pyplot as plt
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from rios import rat
from rios import ratapplier
from tpot import TPOTRegressor
Explanation: Biomass Estimation - Putting the RAT into TPOT
<img src='http://www.auscover.org.au/wp-content/themes/auscover/assets/images/logo.png'>
A demonstration on linking field data to image data using machine learning
DISCLAIMER: This is only a demonstration of some cool data science methods you can pull off on your laptop at home, and not about the science of Biomass estimation. There are some very good scientists and programs looking at how to improve biomass estimation and better quantify the error budget. Whatever comes out of the bottom of this worksheet is to be used for your amusement only :)
_By Peter Scarth (Joint Remote Sensing Research Program) _
Abstract
This is a quick and dirty notebook demonstrating how to link two data sets I pulled off the AusCover portal. In particular, I wanted to show how powerful the combination of Raster Attribute Tables (RATs) and Machine Learning (ML) is for getting quick insights into data at a national scale.  I pulled this together one Sunday night so it's not necessarily pretty or efficient but it might give someone else a heads up into getting started with integrated ecological data science using TERN data.
This uses a number of awesome python libraries. I'd really recommend conda to get going, or if you like Docker pull my datascience-rs container that has everything you need to run this installed.
Raster Attribute Tables
For an excellent introduction to RATs I recommend the following paper:
Clewley, D.; Bunting, P.; Shepherd, J.; Gillingham, S.; Flood, N.; Dymond, J.; Lucas, R.; Armston, J.; Moghaddam, M. A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables. Remote Sensing 2014, 6, 6111-6135. 
RATs work well with machine learning and quick data analysis. To run the prediction across the continent using this notebook takes about 30 minutes (depends on the final model) when working on the RAT. It would take around 6000x longer or nearly 18 weeks to process the image on a per-pixel basis.
If you want a lightweight raster GIS with powerful RAT manipulation abilities have a look at Tuiview
Machine Learning
I really like TPOT for this sort of analysis. It's close to magical for me. You can possibly do better through careful ML pipeline selection and tuning but it really helps to get you started.
"Consider TPOT your Data Science Assistant. TPOT is a Python tool that automatically creates and optimizes machine learning pipelines using genetic programming."
Biomass Library
The Biomass Plot Library is  a national collation of tree and shrub inventory data, allometric model predictions of above and below-ground biomass for Australia. More data is being added all the time.
The data in this notebook was extracted using the Field Data Portal. I selected the Biomass Library and then the Site Level Dataset.  I downloaded the whole lot as a csv file and that file (biolib_sitelist_auscover.csv) is what I use here.
<img src='http://www.auscover.org.au/wp-content/uploads/biomass-plot-library.jpg' >
Height and Structure Data
The Vegetation height and structure - derived from ALOS-1 PALSAR, Landsat and ICESat/GLAS, Australia contains vertical plant profiles for the Australian continent derived through integration of ICESat GLAS waveforms with ALOS PALSAR and Landsat data products. 
There are individual tiff files of the various parameters available for download, but you're going to want the  KEA file of all the parameters and the segments (alpsbk_aust_y2009_sc5a2.kea) which is available here as a 9.1 G download.  I'm yet to work out how to efficiently access kea files over the internet, like we do using the /vsicurl/ driver for GeoTIFF, so for now you'll need to download it.
<img src='http://www.auscover.org.au/wp-content/uploads/Vegetation-height-and-structure.jpg' >
End of explanation
# Read Biomass library data from the csv file 
fieldBiomass=numpy.loadtxt("biolib_sitelist_auscover.csv",delimiter=',', skiprows=1)
# Open Height Map dataset
keaFile = "alpsbk_aust_y2009_sc5a2.kea"
heightDataset = gdal.Open(keaFile, gdal.GA_Update)
# Set up the reprojection transform from WGS84 (biomass library) to Australian Albers (height data)
source = osr.SpatialReference()
source.ImportFromEPSG(4326)
target = osr.SpatialReference()
target.ImportFromEPSG(3577)
transform = osr.CoordinateTransformation(source, target)
# Open the raster band with the segment IDs
heightBand=heightDataset.GetRasterBand(1)
# Get the Albers to pixel transform
geoTransform=heightDataset.GetGeoTransform()
# Find the segmentID for all the field sites
print("Linking field observations to segment IDs\n")
segmentIDs = []
for record in fieldBiomass:
    # Make up a site OGR point
    site = ogr.Geometry(ogr.wkbPoint)
    site.AddPoint(record[0], record[1])
    # Transform the site to EPSG3577
    site.Transform(transform)
    # Get the pixel location of the site
    mx,my=site.GetX(), site.GetY()  #coord in map units
    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - geoTransform[0]) / geoTransform[1]) #x pixel
    py = int((my - geoTransform[3]) / geoTransform[5]) #y pixel
    # Extract the segmentID for the location
    segmentIDs.append(heightBand.ReadAsArray(px,py,1,1)[0][0])
# Get the RAT column names
colNames = rat.getColumnNames(heightDataset)
# Select the columns used for the training/prediction
trainingColumns = [5,6,7,8,9,10,15,19,22,23,25,26,27,28,29,30,31,32,33,34]
trainingNames = [colNames[i] for i in trainingColumns]
# Now we have the segmentIDs, pull the image data from the RAT that corresponds to the segment IDs
imageData = []
# Iterate for all the RAT columns
for name in trainingNames:
    print("Extracting sites from " + name)
    # Extract the array of values corresponding to the field site segments
    imageData.append(rat.readColumnFromBand(heightBand,name).astype('float')[segmentIDs])
# Convert the list of arrays to an array
imageData = numpy.transpose(numpy.array(imageData))
# Remove nodata from the couple of segments too small to get statistics in the image data
goodDataIDX = imageData.min(axis=1)>0
imageData = imageData[goodDataIDX]
fieldBiomass = fieldBiomass[goodDataIDX]
print("\nTraining data has %d observations and %d columns" % imageData.shape)
Explanation: Prepare the training data
This section reads in the Biomass library csv file and links the observations to the closest segment.
NOTE: No consideration is made of the size or shape of the field plot or the shapes of the segments!
The Biomass Library has the following fields:
longitude,latitude,live_basal_area_ha,dead_basal_area_ha,live_tree_density_ha,dead_tree_density_ha,
agb_drymass_ha,agb_drymass_se_ha,bgb_drymass_ha,bgb_drymass_se_ha,tb_drymass_ha,tb_drymass_se_ha  
Once we find the closest segment, we extract the attributes we want from the RAT to aget an imageData array to match with the fieldBiomass array.
End of explanation
# Total amount of time we allow for the training
optTime = 600
# Number of CPUs to use for training and cross validation
nCPUs = 4
# What function to minimise
scoring = 'mean_squared_error'
# Number of subsamples from the Biomass Library for model training
nSubsets = 9999
# We select the tb_drymass_ha column to train on
totalBiomass = fieldBiomass[:,10]
# This is the standard error of the site level estimates
totalBiomassSE = fieldBiomass[:,11]
# Select a subsample to improve the model search speed
subSample = numpy.random.choice(len(totalBiomass),nSubsets,replace=False)
biomass=totalBiomass[subSample]
biomassSE=totalBiomassSE[subSample]
trainData=imageData[subSample]
# Use the proportion of the error in the estimates as fitting weights
biomassWeights=biomass/biomassSE
# Setup the TPOT regression options
tpot = TPOTRegressor(max_time_mins=optTime,  
                     n_jobs = nCPUs, 
                     scoring=scoring, 
                     verbosity=2, 
                     cv=10, 
                     max_eval_time_mins=1,
                     population_size=100)
# Start testing models using 10 fold cross validation and 100 models per generation
tpot.fit(trainData, biomass, sample_weight=biomassWeights)
# Export the best model to a file
tpot.export('tpot_biomass_pipeline.py')
# Build the biomass predictive model
biomassModel = tpot._fitted_pipeline.fit(imageData, totalBiomass)
# Predict the full dataset
predBiomass = biomassModel.predict(imageData)
# Print some RMSE Statistics for various ranges
print("\nTotal RMSE = %f\n" % numpy.sqrt(numpy.mean((totalBiomass-predBiomass)**2)))
stopPoints=[0,100,500,1000,2000,5000,10000]
print("Start"," Stop","Count"," RMSE")
for i in range(len(stopPoints)-1):
    idx=numpy.logical_and(totalBiomass>stopPoints[i],totalBiomass<stopPoints[i+1])
    rmse=numpy.sqrt(numpy.mean((totalBiomass[idx]-predBiomass[idx])**2))
    print('{0:5d} {1:5d} {2:5d} {3:5.0f}'.format(stopPoints[i],stopPoints[i+1],idx.sum(),rmse))
# Plot the Output in a LogLog figure
fig = plt.figure(figsize=(10,10))
plt.loglog(totalBiomass,predBiomass, 'g.',[10,10000], [10,10000],'r-')
plt.xlabel('Observed (Mg/ha)', fontsize=18)
plt.ylabel('Predicted (Mg/ha)', fontsize=18)
plt.title('Total Biomass Estimate', fontsize=32)
plt.xlim([10,10000])
plt.ylim([10,10000])
plt.grid(which='minor', alpha=0.4)                                                
plt.grid(which='major', alpha=0.8)                                                
Explanation: Build the Machine Learning Model
This uses the amazing TPOT tool to find a good machine learning pipeline that predicts biomass. It automates the most tedious part of machine learning by intelligently exploring thousands of possible pipelines to find the best one for your data. There are a few setting below to tweak.
  - optTime is the number of minutes you are prepared to wait for the optimisation. The longer the better.
  - nCPUs is the number of CPUs you have available for the search. You'll also need a minimum of about 2GB of RAM per CPU
  - Scoring is the thing you'll be minimising
  - nSubsets picks a random subset from the training data. Smaller is faster but you might miss some variation in the data
  - The field from the biomass library you want to fit. Here we're selecting the Total Biomass column (number 10)
This will run a genetic algorithm until it times out and selects the best fitting model to use for biomass prediction. It saves the model so you can play with it later if you want.
End of explanation
%%time
def _ratapplier_calc_biomass(info, inputs, outputs):
    
    Calculate Biomass from RAT.
    Called by ratapplier below
    
    ratArray = []
    # Iterate for all the RAT columns
    for name in trainingNames:
        # Extract the array of values corresponding to the field site segments
        ratArray.append(getattr(inputs.inrat, name).astype('float'))
    # Convert the list of arrays to an array
    ratArray = numpy.transpose(numpy.array(ratArray))
    # Predict Biomass
    biomass = biomassModel.predict(ratArray)
    # Make the weird inputs nodata
    biomass[ratArray.min(axis=1) < numpy.finfo(numpy.float32).eps] = 0
    # Save to 'totalBiomass' column (will create if doesn't exist)
    setattr(outputs.outrat,"totalBiomass", biomass)
# Set up ratapplier for input / output
inFile = ratapplier.RatAssociations()
outFile = ratapplier.RatAssociations()
# Pass in clumps file for the input and output as we'll be updating the existing RAT
inFile.inrat = ratapplier.RatHandle(keaFile)
outFile.outrat = ratapplier.RatHandle(keaFile)
# Apply function to all rows in chunks
ratapplier.apply(_ratapplier_calc_biomass, inFile, outFile)
    
    
Explanation: Predict Biomass using the RAT
Once we have the biomassModel we can use that to predict the Biomass for every entry in the RAT. 
This is cool, because for Australia we only need to make around 3,000,000 predictions which is far faster than predicting every one of the 18,407,102,543 pixels in the image. That's around a 6000x speedup.
Here,  RIOS ratapplier does the work writing a brand new attribute column into the kea file.
End of explanation
# Additional Imports
import rsgislib
from rsgislib.rastergis import exportCol2GDALImage
# Setup the export parameters
outimage='/home/jovyan/work/Temp/totalBiomass.tif'
gdalformat = 'GTIFF'
datatype = rsgislib.TYPE_16UINT
field = 'totalBiomass'
# Run the export
exportCol2GDALImage(keaFile, outimage, gdalformat, datatype, field)
Explanation: Export the Biomass band to an individual GeoTIFF
You can run this section if you want to output just the newly calculated biomass column to a separate GeoTIFF. 
This uses one of the function from the rastergis module in RSGISLIB
After this, you might like to run gdal_translate to compress the output and gdaladdo to add overviews for easy storage and fast viewing in your favorite GIS platform. Like this:
```
gdal_translate -a_nodata 0 -of GTIFF -co COMPRESS=DEFLATE -co ZLEVEL=9 -co TILED=YES totalBiomass.tif totalBiomassCompressed.tif
gdaladdo -r average --config TILED_OVERVIEW YES --config COMPRESS_OVERVIEW PACKBITS totalBiomassCompressed.tif 4 8 16 32 64 128 256 512
```
If everything has worked, you should have a file that looks something like this:
<img src='https://pbs.twimg.com/media/C-s2iyTUwAAxNzb.jpg' >
End of explanation
<END_TASK>