interpol_comp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Jupyter Notebook is intended to reduce the amount of time spent manually comparing results from interpolation tools in ArcGIS Pro. The notebook was created in May 2020. Questions can be sent to Kalle Jahn at [kjahn@contractor.usgs.gov](mailto:kjahn@contractor.usgs.gov)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import arcpy\n",
    "import itertools\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xml.etree.ElementTree as ET\n",
    "from datetime import datetime\n",
    "\n",
    "arcpy.env.overwriteOutput=True\n",
    "\n",
    "# To enable the extension \"GeoStats\", this seems required for concurrent license users, \n",
    "# e.g., some university users\n",
    "arcpy.CheckOutExtension(\"GeoStats\")\n",
    "arcpy.CheckOutExtension(\"Spatial\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# geodatabase within the project folder that contains your point feature classes\n",
    "projGeodb = 'test4.gdb'\n",
    "arcpy.env.workspace = projGeodb\n",
    "\n",
    "# point feature class name\n",
    "inPoints = 'Combined_L1_L2'\n",
    "\n",
    "# percentage training points\n",
    "trainSize = 90\n",
    "\n",
    "# data to be interpolated (column in attribute table)\n",
    "zField = 'L1_Surface_Meters'\n",
    "\n",
    "# select transformations for simple, ordinary, and universal kriging\n",
    "# TransformationNormalScore is for simple only...\n",
    "# you must at the very least have 'NONE' ...\n",
    "krigTransform = ['TransformationNormalScore','NONE','TransformationBoxCox','TransformationLog']\n",
    "\n",
    "# EBK: neighborhood radius taken from Geostat Wizard\n",
    "ebkRadius = 205\n",
    "\n",
    "# EBK: maximum number of points in local models\n",
    "# note that larger max leads to longer runtime, but necessary for large datasets (>10,000 points)\n",
    "maxLocalPoints = 100\n",
    "\n",
    "# set snap raster file, or write None if you have no snap raster\n",
    "snapRaster = None\n",
    "arcpy.env.snapRaster = snapRaster\n",
    "\n",
    "# if snap raster exits, raster cell size is pulled from snap raster\n",
    "# otherwise generated rasters will have cell size entered manually below\n",
    "manualCellSize = 38.1 # # MUST BE SAME UNITS AS THOSE IN POINT FEATURE CLASS!\n",
    "    \n",
    "# boundary for topo to raster\n",
    "# set to file name, otherwise set to None if you aren't using a boundary\n",
    "boundary = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pre-transformation data shift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shift_up(in_features,z_field):\n",
    "    data = arcpy.da.FeatureClassToNumPyArray(in_features, field_names=[z_field],\n",
    "                                                     skip_nulls=True)\n",
    "    minimum = np.min(data[zField])\n",
    "\n",
    "    if minimum <= 0.0:\n",
    "        factor = 0.1 if minimum == 0 else abs(minimum)+0.1\n",
    "        # for some reasons, I can't get the original code working without adding the column first \n",
    "        # before calculating the field\n",
    "        # My ArcGIS Pro version is 2.4\n",
    "        arcpy.AddField_management(in_features, 'shifted','DOUBLE')\n",
    "        arcpy.CalculateField_management(in_features, 'shifted', f'!{z_field}! + {factor}')\n",
    "    \n",
    "    else:\n",
    "        factor = None\n",
    "    return factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Updating XML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def updateKrigXML(filename,transform,semivariogram,optimize='ByCrossvalidation'):\n",
    "    '''\n",
    "    Updates the XML file with the given semivariogram and adds optimization as if\n",
    "    user had pressed the \"Optimize\" button in Geostat Wizard.\n",
    "    \n",
    "    INPUTS:\n",
    "     filename (XML file):\n",
    "         The XML file to be updated.\n",
    "     semivariogram:\n",
    "         The desired semivariogram model type. Must be available in Geostat Wizard...\n",
    "    OUTPUTS:\n",
    "     Updated XML file.\n",
    "    '''\n",
    "    tree = ET.ElementTree(file=filename)\n",
    "    root = tree.getroot()\n",
    "    \n",
    "    # add the optimize attribute\n",
    "    root.set('optimize', optimize)\n",
    "    \n",
    "    # get krig type\n",
    "    km = root.find('enum[@name=\"KrigingMethodType\"]')\n",
    "    dataset = tree.find('./items/item')\n",
    "\n",
    "    # edit transform type\n",
    "    tran = dataset.find('./model')\n",
    "    tran.set('name',transform)\n",
    "    param = tran.find('./value')\n",
    "    if tran.attrib['name'] == 'TransformationBoxCox' and param == None:\n",
    "        p = ET.SubElement(tran, \"value\",attrib={'name':'Parameter'}).text = '1'  \n",
    "    \n",
    "    # if simple kriging and None or Log or BoxCox, set the mean to auto\n",
    "#     if km.text == 'Simple' and (tran.attrib['name'] == 'None' or tran.attrib['name'] == 'TransformationBoxCox' or tran.attrib['name'] == 'TransformationLog'):\n",
    "    if km.text == 'Simple' and tran.attrib['name'] != 'TransformationNormalScore':\n",
    "        mean = dataset.find('./value[@name=\"Mean\"]')\n",
    "        if mean == None:\n",
    "            ET.SubElement(dataset, \"value\",attrib={'name':'Mean','auto':'true'}).text = '2.0'\n",
    "        else:\n",
    "            mean.set('auto','true')\n",
    "                          \n",
    "    # add \"auto\" to the DensitySkew if that is what was chosen.\n",
    "    t = tree.find('./items/item/model')\n",
    "    m = t.find('./enum[@name=\"ApproximationMethodType\"]')\n",
    "    if m != None and m.text == 'DensitySkew':\n",
    "        t.find('./enum[@name=\"BaseDistribution\"]').set('auto','true')\n",
    "        t.find('./value[@name=\"Kernels\"]').set('auto','true')\n",
    "    \n",
    "    # set the variogram model\n",
    "    sm = root.find('./model/model/enum[@name=\"ModelType\"]')\n",
    "    sm.text = semivariogram\n",
    "    \n",
    "    # need to add another parameter for some of the models\n",
    "    vm = root.find('./model/model[@name=\"VariogramModel\"]')\n",
    "    param = root.find('./model/model/value[@name=\"Parameter\"]')\n",
    "    atts = {\"auto\": \"true\", \"name\": \"Parameter\"}\n",
    "    if semivariogram == 'Stable':\n",
    "        if param == None:\n",
    "            new = ET.SubElement(vm, \"value\",attrib=atts).text = '2'  \n",
    "        else:\n",
    "            param.text = '2'\n",
    "    elif semivariogram == 'BesselK' or semivariogram == 'BesselJ':\n",
    "        if param == None:\n",
    "            new = ET.SubElement(vm, \"value\",attrib=atts).text = '9.999'\n",
    "        else:\n",
    "            param.text = '9.999'\n",
    "    else:\n",
    "        if param != None:\n",
    "            vm.remove(param)\n",
    "        \n",
    "    \n",
    "    with open(filename, 'wb') as fh:\n",
    "        tree.write(fh)\n",
    "    fh.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross validation of Kriging and EBK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crossVal(outLayer): \n",
    "    '''\n",
    "    INPUTS:\n",
    "     outLayer (Geostatistical Layer):\n",
    "         The geostatistical layer to be analyzed.\n",
    "    OUTPUTS:\n",
    "     statistics list :\n",
    "         The cross-validation statistics stored in list form.\n",
    "    '''\n",
    "    cvResult = arcpy.CrossValidation_ga(outLayer)\n",
    "    \n",
    "    meanErr = float(cvResult.meanError)\n",
    "    meanStd = float(cvResult.meanStandardized)\n",
    "    rmsErr = float(cvResult.rootMeanSquare)\n",
    "    rmsStd = float(cvResult.rootMeanSquareStandardized)\n",
    "    avgStd = float(cvResult.averageStandard)\n",
    "    in90 = float(cvResult.percentIn90Interval)\n",
    "    in95 = float(cvResult.percentIn95Interval)\n",
    "    avgCRPS = float(cvResult.averageCRPS)\n",
    "    \n",
    "    return [outLayer,meanErr,meanStd,rmsErr,rmsStd,avgStd,in90,in95,avgCRPS]        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subsetting, Training, and Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subset(in_features,trainSize=90):\n",
    "    '''\n",
    "    Subset the input points to allow for training and testing comparison.\n",
    "    \n",
    "    INPUTS:\n",
    "     in_features (points):\n",
    "         Point to be kriged.\n",
    "     trainSize:\n",
    "         Percentage of points to use in the training set.\n",
    "    OUTPUTS:\n",
    "     trainPoints:\n",
    "        Points for creating the interpolated surface.\n",
    "     testPoints:\n",
    "        Points for testing the interpolated surface.\n",
    "    \n",
    "    '''\n",
    "    trainPoints = f'{in_features[:-4]}_train_pts{trainSize}'\n",
    "    testPoints = f'{in_features[:-4]}_test_pts{100-trainSize}'\n",
    "    subsizeUnits = \"PERCENTAGE_OF_INPUT\"\n",
    "\n",
    "    # check to see if the trainPoints and testPoints already exist\n",
    "    # delete them if they do\n",
    "    for fc in [trainPoints, testPoints]:\n",
    "        if arcpy.Exists(fc):\n",
    "            arcpy.Delete_management(fc)\n",
    "            \n",
    "    arcpy.SubsetFeatures_ga(in_features, \n",
    "                            trainPoints, \n",
    "                            testPoints,\n",
    "                            trainSize, \n",
    "                            subsizeUnits)\n",
    "    return trainPoints, testPoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TrainTestKrig(train_points,\n",
    "                  test_points,\n",
    "                  z_field,\n",
    "                  outLayer,\n",
    "                  xml,\n",
    "                  cell_size,\n",
    "                  trainSize=90\n",
    "                 ):\n",
    "    '''\n",
    "    TO DO\n",
    "    \n",
    "    '''\n",
    "    # create GA layers with 90% training points\n",
    "    # method for Kriging\n",
    "    outLayerTrain = f'{outLayer}_train'\n",
    "    newDatasetTrain = arcpy.GeostatisticalDatasets(xml)\n",
    "    newDatasetTrain.dataset1 = train_points\n",
    "    newDatasetTrain.dataset1Field = z_field\n",
    "    arcpy.GACreateGeostatisticalLayer_ga(xml,newDatasetTrain, outLayerTrain)\n",
    "    \n",
    "#     # save outLayer to GA layer file\n",
    "#     out_layer_file = '.\\\\geostat_layers\\\\'+outLayerTrain+'.lyrx'\n",
    "#     arcpy.SaveToLayerFile_management(outLayerTrain, out_layer_file, 'RELATIVE')\n",
    "\n",
    "    # GA to raster to points using outLayerTrain and 90% train points\n",
    "    trainRaster = f'{outLayer}_train_raster'\n",
    "    try:\n",
    "        arcpy.GALayerToGrid_ga(outLayerTrain, trainRaster, cell_size=cell_size)\n",
    "        # extract raster values\n",
    "        trainResult = f'{outLayer}_train_result'\n",
    "        arcpy.sa.ExtractValuesToPoints(train_points, trainRaster, trainResult)\n",
    "        # export values to numpy array\n",
    "        field_names = [z_field,'RASTERVALU']\n",
    "        arrTrain = arcpy.da.FeatureClassToNumPyArray(trainResult, field_names=field_names,\n",
    "                                                         skip_nulls=False)\n",
    "        meanTrain = np.mean(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "        medianTrain = np.median(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "        arcpy.Delete_management(trainResult)\n",
    "    except arcpy.ExecuteError:\n",
    "        meanTrain = f'ExecuteError {arcpy.getMessages()}'\n",
    "        medianTrain = f'ExecuteError {arcpy.getMessages()}'\n",
    "    \n",
    "    # GA to points using trainRaster and 10% test points\n",
    "    try:\n",
    "        # extract raster values\n",
    "        testResult = f'{outLayer}_test_result'\n",
    "        x = 'extract values'\n",
    "        arcpy.sa.ExtractValuesToPoints(test_points, trainRaster, testResult)\n",
    "        # export values to numpy array\n",
    "        field_names = [z_field,'RASTERVALU']\n",
    "        x = 'to array'\n",
    "        arrTest = arcpy.da.FeatureClassToNumPyArray(testResult, field_names=field_names,\n",
    "                                                         skip_nulls=False)\n",
    "        meanTest = np.mean(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "        medianTest = np.median(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "        arcpy.Delete_management(testResult)\n",
    "    except arcpy.ExecuteError:\n",
    "        meanTest = f'ExecuteError {arcpy.getMessages()}'\n",
    "        medianTest = f'ExecuteError {arcpy.getMessages()}'\n",
    "    \n",
    "    arcpy.Delete_management(outLayerTrain)\n",
    "    if trainRaster != None: arcpy.Delete_management(trainRaster)\n",
    "        \n",
    "    # return results in list\n",
    "    return [outLayer,meanTrain,medianTrain,meanTest,medianTest]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TrainTestEBK(train_points,\n",
    "                 test_points,\n",
    "                 z_field,\n",
    "                 outLayer,\n",
    "                 outLayerTrain,\n",
    "                 cell_size\n",
    "                ):\n",
    "    '''\n",
    "    TO DO\n",
    "    \n",
    "    '''    \n",
    "    # GA to raster and points using outLayerTrain and 90% train points\n",
    "    trainRaster = f'{outLayer}_train_raster'    \n",
    "    try:\n",
    "        arcpy.GALayerToGrid_ga(outLayerTrain,trainRaster,cell_size=cell_size)\n",
    "        # extract raster values\n",
    "        trainResult = f'{outLayer}_train_result'\n",
    "        arcpy.sa.ExtractValuesToPoints(train_points,trainRaster,trainResult)\n",
    "        # export values to numpy array\n",
    "        field_names = [z_field,'RASTERVALU']\n",
    "        arrTrain = arcpy.da.FeatureClassToNumPyArray(trainResult, field_names=field_names,\n",
    "                                                         skip_nulls=False)\n",
    "        meanTrain = np.mean(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "        medianTrain = np.median(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "        arcpy.Delete_management(trainResult)\n",
    "    except arcpy.ExecuteError:\n",
    "        meanTrain = f'ExecuteError {arcpy.getMessages()}'\n",
    "        medianTrain = f'ExecuteError {arcpy.getMessages()}'\n",
    "\n",
    "    # GA to raster to points using trainRaster and 10% test points\n",
    "    try:\n",
    "        # extract raster values\n",
    "        testResult = f'{outLayer}_test_result'\n",
    "        x = 'extract values'\n",
    "        arcpy.sa.ExtractValuesToPoints(test_points,trainRaster,testResult)\n",
    "        # export values to numpy array\n",
    "        field_names = [z_field,'RASTERVALU']\n",
    "        x = 'to array'\n",
    "        arrTest = arcpy.da.FeatureClassToNumPyArray(testResult, field_names=field_names,\n",
    "                                                         skip_nulls=False)\n",
    "        meanTest = np.mean(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "        medianTest = np.median(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "        arcpy.Delete_management(testResult)\n",
    "    except arcpy.ExecuteError:\n",
    "        meanTest = f'ExecuteError {arcpy.getMessages()}'\n",
    "        medianTest = f'ExecuteError {arcpy.getMessages()}'\n",
    "    if trainRaster != None: arcpy.Delete_management(trainRaster)\n",
    "    # return results in list\n",
    "    return [outLayer,meanTrain,medianTrain,meanTest,medianTest]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automated Kriging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autoKrig(in_features,\n",
    "             train_points,\n",
    "             test_points,\n",
    "             z_field,\n",
    "             cell_size,\n",
    "             transform = ['TransformationNormalScore','NONE','TransformationBoxCox','TransformationLog'],\n",
    "             krigTypes = ['Simple','Ordinary','Universal'],\n",
    "             semivariograms = ['Circular','Spherical','Exponential',\n",
    "                               'Gaussian','Stable','BesselK',\n",
    "                               'BesselJ','RationalQuadratic'],\n",
    "             ):\n",
    "    '''\n",
    "    TO ADD \n",
    "    \n",
    "    '''\n",
    "    # empty lists to fill with results\n",
    "    k_results = []\n",
    "    subset_results = []\n",
    "    \n",
    "    # models broken up below because only simple kriging has NormalScore transform\n",
    "    simModels = list(itertools.product(['Simple'],transform,semivariograms))\n",
    "    if 'TransformationNormalScore' in transform: \n",
    "        new = transform.copy()\n",
    "        new.remove('TransformationNormalScore')\n",
    "    else:\n",
    "        new = transform.copy()\n",
    "    ordUnvModels = list(itertools.product(['Ordinary','Universal'],new,semivariograms))\n",
    "    \n",
    "    models = simModels + ordUnvModels\n",
    "    \n",
    "    for i in models:\n",
    "        if i[0] == 'Ordinary':\n",
    "            xml = 'KrigingOrd.xml'\n",
    "        elif i[0] == 'Simple':\n",
    "            xml = 'KrigingSim.xml'\n",
    "        elif i[0] == 'Universal':\n",
    "            xml = 'KrigingUnv.xml'\n",
    "            \n",
    "        outLayer = f'Kriging_{i[0]}_{i[1]}_{i[2]}'\n",
    "        print(f'{datetime.now().strftime('%H:%M:%S')}: working on {outLayer}')\n",
    "\n",
    "        # update the XML file with appropriate semivariogram and transformation \n",
    "        updateKrigXML(xml,i[1],i[2])\n",
    "\n",
    "        # create new GA dataset\n",
    "        newDataset = arcpy.GeostatisticalDatasets(xml)\n",
    "        newDataset.dataset1 = in_features\n",
    "        newDataset.dataset1Field = z_field\n",
    "        # new GA layer with optimized parameters\n",
    "        arcpy.GACreateGeostatisticalLayer_ga(xml, newDataset, outLayer)\n",
    "\n",
    "        # save outLayer to GA layer file\n",
    "        out_layer_file = f'./{inPoints[:-4]}_geostat_layers/{outLayer}.lyrx'\n",
    "        arcpy.SaveToLayerFile_management(outLayer, out_layer_file, 'RELATIVE')\n",
    "\n",
    "        # calculate statistics and append to list\n",
    "        try:\n",
    "            r = crossVal(outLayer)\n",
    "            k_results.append(r)\n",
    "        except:\n",
    "            k_results.append(['CrossValExeErr'] * 9)\n",
    "        # train and test with subsets            \n",
    "        sr = TrainTestKrig(train_points,test_points,z_field,outLayer,xml,cell_size)\n",
    "        subset_results.append(sr)\n",
    "\n",
    "        arcpy.Delete_management(outLayer)\n",
    "    \n",
    "    print('creating dataframes')\n",
    "    # create pandas dataframes from lists\n",
    "    df1 = pd.DataFrame(k_results,columns=['model','ME','MSDE','RMSE','RMSSDE','avgStd','in90','in95','avgCRPS'])\n",
    "    df2 = pd.DataFrame(subset_results,columns=['model','Mean Err - training pts','Median Err - training pts',\n",
    "                                                       'Mean Err - testing pts','Median Err - testing pts'])\n",
    "    return df1, df2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automated EBK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# EBK initial test\n",
    "def autoEBK(in_features,\n",
    "            train_points,\n",
    "            test_points,\n",
    "            z_field,\n",
    "            cell_size,\n",
    "            ebk_radius,\n",
    "            maxLocalPoints,\n",
    "            numberSemivariograms=100\n",
    "           ):\n",
    "    '''\n",
    "    TO ADD \n",
    "    \n",
    "    '''\n",
    "    \n",
    "    field_names = [z_field]\n",
    "    checkNeg = arcpy.da.FeatureClassToNumPyArray(in_features, field_names=field_names,\n",
    "                                                     skip_nulls=True)\n",
    "    minimum = np.min(checkNeg[zField])\n",
    "    transformations = ['NONE','EMPIRICAL'] if minimum <=0 else ['NONE', 'EMPIRICAL', 'LOGEMPIRICAL']\n",
    "#     transformations = ['NONE']\n",
    "    \n",
    "    # Search neighborhood variables\n",
    "    # Standard Circular\n",
    "    angle = 0\n",
    "    maxNeighbors = 15\n",
    "    minNeighbors = 5\n",
    "    sectorType = 'FOUR_SECTORS_SHIFTED'\n",
    "    searchNeighbourhood = arcpy.SearchNeighborhoodStandardCircular(ebk_radius, angle, maxNeighbors,\n",
    "                                                                   minNeighbors, sectorType)\n",
    "    # parameters\n",
    "    outputType = 'PREDICTION'\n",
    "    quantileValue = None\n",
    "    thresholdType = None\n",
    "    probabilityThreshold = None\n",
    "    overlapFactor = 1\n",
    "    \n",
    "    ebk_results = []\n",
    "    subset_results = []\n",
    "    \n",
    "    # these semivariograms are for no transformation\n",
    "    s1 = ['POWER','LINEAR','THIN_PLATE_SPLINE']\n",
    "    # these semivariograms are for empirical and logempirical transformations\n",
    "    s2 = ['EXPONENTIAL','EXPONENTIAL_DETRENDED','WHITTLE','WHITTLE_DETRENDED','K_BESSEL','K_BESSEL_DETRENDED']\n",
    "    \n",
    "    # create model combinations list\n",
    "    models = list(itertools.product(transformations[0:1],s1)) + \\\n",
    "             list(itertools.product(transformations[1:],s2))\n",
    "    \n",
    "    for i in models:\n",
    "        outLayer = f'EBK_{i[0]}_{i[1]}'\n",
    "        print(f'{datetime.now().strftime('%H:%M:%S')}: working on {outLayer}')\n",
    "        \n",
    "        arcpy.EmpiricalBayesianKriging_ga(in_features=in_features,\n",
    "                                  z_field=z_field,\n",
    "                                  out_ga_layer=outLayer,\n",
    "                                  out_raster=None,\n",
    "                                  cell_size=None,\n",
    "                                  transformation_type=i[0],\n",
    "                                  max_local_points=maxLocalPoints,\n",
    "                                  overlap_factor=overlapFactor,\n",
    "                                  number_semivariograms=numberSemivariograms,\n",
    "                                  search_neighborhood=searchNeighbourhood, \n",
    "                                  output_type=outputType, \n",
    "                                  quantile_value=quantileValue, \n",
    "                                  threshold_type=thresholdType, \n",
    "                                  probability_threshold=probabilityThreshold,\n",
    "                                  semivariogram_model_type=i[1]\n",
    "                                 )\n",
    "        \n",
    "        # save GA layer file\n",
    "        out_layer_file = f'./{inPoints[:-4]}_geostat_layers/{outLayer}.lyrx'\n",
    "        arcpy.SaveToLayerFile_management(outLayer, out_layer_file, 'RELATIVE')\n",
    "\n",
    "        # Cross-validation\n",
    "        r = crossVal(outLayer)\n",
    "        ebk_results.append(r)\n",
    "\n",
    "        # train and test with subsets\n",
    "        # create training layer\n",
    "        outLayerTrain = f'{outLayer}_train'\n",
    "        arcpy.EmpiricalBayesianKriging_ga(in_features=train_points,\n",
    "                                          z_field=z_field,\n",
    "                                          out_ga_layer=outLayerTrain,\n",
    "                                          out_raster=None,\n",
    "                                          cell_size=None,\n",
    "                                          transformation_type=i[0],\n",
    "                                          max_local_points=maxLocalPoints,\n",
    "                                          overlap_factor=overlapFactor,\n",
    "                                          number_semivariograms=numberSemivariograms,\n",
    "                                          search_neighborhood=searchNeighbourhood, \n",
    "                                          output_type=outputType, \n",
    "                                          quantile_value=quantileValue, \n",
    "                                          threshold_type=thresholdType, \n",
    "                                          probability_threshold=probabilityThreshold,\n",
    "                                          semivariogram_model_type=i[1]\n",
    "                                         )\n",
    "\n",
    "        sr = TrainTestEBK(train_points,test_points,z_field,outLayer,outLayerTrain,cell_size)\n",
    "        subset_results.append(sr)\n",
    "        \n",
    "        arcpy.Delete_management(outLayer)\n",
    "        arcpy.Delete_management(outLayerTrain)\n",
    "                \n",
    "    print('creating dataframes')\n",
    "    # create pandas dataframes from lists\n",
    "    df1 = pd.DataFrame(ebk_results,columns=['model','ME','MSDE','RMSE','RMSSDE','avgStd','in90','in95','avgCRPS'])\n",
    "    df2 = pd.DataFrame(subset_results,columns=['model','Mean Err - training pts','Median Err - training pts',\n",
    "                                                       'Mean Err - testing pts','Median Err - testing pts'])\n",
    "\n",
    "    \n",
    "    return df1, df2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Topo to raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ttr(trainPoints,\n",
    "        testPoints,\n",
    "        z_field,\n",
    "        cell_size,\n",
    "        boundary\n",
    "       ):\n",
    "        \n",
    "    inPointElevations = f'{trainPoints} {z_field} POINTELEVATION'\n",
    "    if boundary == None:\n",
    "        outRaster = f'{inPoints[:-4]}_topo_no_bound'\n",
    "        arcpy.TopoToRaster_3d(inPointElevations, outRaster, cell_size=cell_size, data_type='SPOT')\n",
    "    else:\n",
    "        outRaster = inPoints[:-4]+f'{inPoints[:-4]}_topo_bound'\n",
    "        inBoundary = f'{boundary} # BOUNDARY'\n",
    "        inFeats = (f'{inPointElevations};{inBoundary}')\n",
    "        arcpy.TopoToRaster_3d(inFeats, outRaster, cell_size=cell_size, data_type='SPOT') \n",
    "\n",
    "    # compare raster to 90% test points\n",
    "    trainResult = f'{inPoints[:-4]}_topoTrain'\n",
    "    arcpy.sa.ExtractValuesToPoints(trainPoints,outRaster,trainResult)\n",
    "    \n",
    "    field_names = [z_field,'RASTERVALU']\n",
    "    arrTrain = arcpy.da.FeatureClassToNumPyArray(trainResult, field_names=field_names,\n",
    "                                                     skip_nulls=True)\n",
    "    meanTrain = np.mean(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "    medianTrain = np.median(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "    \n",
    "    # compare raster to 10% test points\n",
    "    testResult = f'{inPoints[:-4]}_topoTest'\n",
    "    arcpy.sa.ExtractValuesToPoints(testPoints,outRaster,testResult)\n",
    "    \n",
    "    arrTest = arcpy.da.FeatureClassToNumPyArray(testResult, field_names=field_names,\n",
    "                                                     skip_nulls=True)\n",
    "    meanTest = np.mean(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "    medianTest = np.median(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "    \n",
    "    arcpy.Delete_management(outRaster)\n",
    "    arcpy.Delete_management(trainResult)\n",
    "    arcpy.Delete_management(testResult)\n",
    "    \n",
    "    return [outRaster,meanTrain,medianTrain,meanTest,medianTest]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automated raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autoRaster(trainPoints,\n",
    "               testPoints,\n",
    "               z_field,\n",
    "               cell_size,\n",
    "               boundary\n",
    "              ):\n",
    "        \n",
    "    results = []\n",
    "    \n",
    "    # NATURAL NEIGHBORS\n",
    "    print(datetime.now().strftime('%H:%M:%S')+' working on Natural Neighbor')\n",
    "    outNat = arcpy.sa.NaturalNeighbor(trainPoints,z_field,cell_size)\n",
    "    \n",
    "    # compare raster to 90% test points\n",
    "    trainResult = inPoints[:-4]+'_natTrain'\n",
    "    arcpy.sa.ExtractValuesToPoints(trainPoints,outNat,trainResult)\n",
    "    \n",
    "    field_names = [z_field,'RASTERVALU']\n",
    "    arrTrain = arcpy.da.FeatureClassToNumPyArray(trainResult, field_names=field_names,\n",
    "                                                     skip_nulls=True)\n",
    "    meanTrain = np.mean(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "    medianTrain = np.median(abs(arrTrain[z_field] - arrTrain['RASTERVALU']))\n",
    "    \n",
    "    # compare raster to 10% test points\n",
    "    testResult = inPoints[:-4]+'_natTest'\n",
    "    arcpy.sa.ExtractValuesToPoints(testPoints,outNat,testResult)\n",
    "    \n",
    "    arrTest = arcpy.da.FeatureClassToNumPyArray(testResult, field_names=field_names,\n",
    "                                                     skip_nulls=True)\n",
    "    meanTest = np.mean(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "    medianTest = np.median(abs(arrTest[z_field] - arrTest['RASTERVALU']))\n",
    "\n",
    "    results.append(['Natural_Neighbor',meanTrain,medianTrain,meanTest,medianTest])\n",
    "    \n",
    "    arcpy.Delete_management(trainResult)\n",
    "    arcpy.Delete_management(testResult)\n",
    "    \n",
    "    # TOPO TO RASTER\n",
    "    print(datetime.now().strftime('%H:%M:%S')+' working on Topo to Raster without boundary')\n",
    "    nobound = ttr(trainPoints,testPoints,z_field,cell_size,boundary=None)\n",
    "    results.append(nobound)\n",
    "    \n",
    "    if boundary != None:\n",
    "        print(datetime.now().strftime('%H:%M:%S')+' working on Topo to Raster with boundary')\n",
    "        bound = ttr(trainPoints,testPoints,z_field,cell_size,boundary)\n",
    "        results.append(bound)\n",
    "        \n",
    "    # prepare dataframe\n",
    "    df = pd.DataFrame(results,columns=['model','Mean Err - training pts','Median Err - training pts',\n",
    "                                          'Mean Err - testing pts','Median Err - testing pts'])\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Execution cells"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data shift, train/test points, set cell size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalize the data set if necessary and change zField if shift occurs\n",
    "factor = shift_up(inPoints,zField)\n",
    "print('normalize factor = {}'.format(factor))\n",
    "if factor != None:\n",
    "    zField = 'shifted'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create training and testing points\n",
    "trainPoints, testPoints = subset(inPoints, trainSize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select cell size based on snapraster\n",
    "if snapRaster != None:\n",
    "    cellSize = float(arcpy.GetRasterProperties_management(snapRaster, 'CELLSIZEY').getOutput(0))\n",
    "else:\n",
    "    cellSize = manualCellSize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Kriging methods (remove Bessel models if they create an error)\n",
    "k1, k2 = autoKrig(inPoints,trainPoints,testPoints,zField,cellSize,krigTransform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Empirical Bayesian Kriging\n",
    "ebk1, ebk2 = autoEBK(inPoints,trainPoints,testPoints,zField,cellSize,ebkRadius,maxLocalPoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Topo to Raster and Natural Neighbors\n",
    "raster = autoRaster(trainPoints,testPoints,zField,cellSize,boundary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DATAFRAME CONCATINATION \n",
    "stats = pd.concat([k1,ebk1],sort=False)#.sort_values(by='RMSE')\n",
    "sub_tests = pd.concat([k2,ebk2,raster],sort=False)#.sort_values(by='Mean Err - testing pts')\n",
    "f = pd.DataFrame([factor], columns=['factor'],index=[''])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DATAFRAME EXPORT TO EXCEL\n",
    "with pd.ExcelWriter(f'{inPoints[:-4]}_interpol.xlsx') as writer: # use this to create new file\n",
    "    stats.to_excel(writer, sheet_name='Statistics')\n",
    "    sub_tests.to_excel(writer, sheet_name='Subset_comparison')\n",
    "    f.to_excel(writer, sheet_name='normalize_factor')\n",
    "writer.save()\n",
    "print('saved to Excel file')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for concurrent users, return the extensions\n",
    "arcpy.CheckInExtension(\"GeoStats\")\n",
    "arcpy.CheckInExtension(\"Spatial\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}