PHP Code

The PHP program used to assign treatment prescriptions is broken into two files:

  • AssignTreatments.php
  • CentroidMatrix.php
The source code for AssignTreatments.php:
<?php
// Program to read a CSV file of tamarisk presence polygons, and
// assign treatments to each polygon based on attributes
//
// Ted Manahan
// NR 505 Fall 2010
require("./CentroidMatrix.php");

// The lowest average polygon cost (best accessibility).
// This is used to scale treatment costs relative to accessibility
// Future improvement: calculate this from input data

// For the Lower Gunnison River
define("BASECOST",148013.7031);

// For the John Martin SWA
//define("BASECOST",1005);

//error_reporting(E_ALL);
ini_set('display_errors', '1');
date_default_timezone_set("America/Denver");

// Some global threshold numbers to use when assigning treatments
// Aerial control is only appropriate for dense stands
$MinPctCovForAerial=75;
// Aerial control is only appropriate for large stands
// 250 acres expressed in square meters
$MinAreaForAerial=250*4046.85642;
// Mechanical only works on fairly flat ground
$MaxSlopeForMechanical=20;
// Mechanical only worthwhile for somewhat dense stands
$MinPctCovForMechanical=50;

// function to aggregate treatment groups
function AggregateTGs ($TreatmentGroups, $TamariskPolygons)
{
    // Default values to use for missing data
    $DefaultHeight=2;
    $DefaultPct_Cov=50;

    // The number of treatment groups
    $TGCount=count($TreatmentGroups);

    // Aggregate treatment groups
    $Group=array();
    for ($TG=0; $TG<$TGCount; $TG++)
    {
        $Group[$TG]["Area"]=0;
        $Group[$TG]["SlopeMAX"]=0;
        $Group[$TG]["Height"]=0;
        $Group[$TG]["Pct_Cov"]=0;
        $Group[$TG]["CostMEAN"]=0;

        $NumPolysInTG=count($TreatmentGroups[$TG]);
        for ($TGPoly=0;$TGPoly<$NumPolysInTG;$TGPoly++)
        {
            $Poly=$TreatmentGroups[$TG][$TGPoly];

            // Get total area of all polygons in group
            $Group[$TG]["Area"]+=$TamariskPolygons[$Poly]["Area"];

            // Look for the maximum slope over all the polygons in group
            if ($TamariskPolygons[$Poly]["SlopeMAX"]>
                    $Group[$TG]["SlopeMAX"])
                $Group[$TG]["SlopeMAX"]=$TamariskPolygons[$Poly]["SlopeMAX"];

            // Sum height, calculate average later
            $Height=$TamariskPolygons[$Poly]["Height"];
            if ($Height<=0)$Height=$DefaultHeight;
            $Group[$TG]["Height"]+=$Height;

            // Sum cover percents, calculate average later
            $Pct_Cov=$TamariskPolygons[$Poly]["Pct_Cov"];
            if ($Pct_Cov<=0)$Pct_Cov=$DefaultPct_Cov;
            $Group[$TG]["Pct_Cov"]+=$Pct_Cov;

            // Sum costs, calculate mean later
            $Group[$TG]["CostMEAN"]+=$TamariskPolygons[$Poly]["CostMEAN"];
        }
       
        // Averages over treatment group
        $Group[$TG]["Height"]  =$Group[$TG]["Height"]  / $NumPolysInTG;
        $Group[$TG]["Pct_Cov"] =$Group[$TG]["Pct_Cov"] / $NumPolysInTG;
        $Group[$TG]["CostMEAN"]=$Group[$TG]["CostMEAN"]/ $NumPolysInTG;
    }
    return($Group);
}

// function to assign the treatments
function AssignTreatment($TamariskPolygons, $TreatmentGroups)
{
    global $MinPctCovForAerial;
    global $MinAreaForAerial;
    global $MaxSlopeForMechanical;
    global $MinPctCovForMechanical;

    // The number of polygons and treatment groups
    $NumPolygons=count($TamariskPolygons);
    $TGCount=count($TreatmentGroups);
    echo "AssignTreatment: Assigning treatments for $NumPolygons polygons ".
         "in $TGCount groups.\n";

    // Aggregate treatment groups
    //echo "AssignTreatment: Calling AggregateTGs from AssignTreatment\n";
    $GroupAttrs=AggregateTGs ($TreatmentGroups, $TamariskPolygons);

    // Assign treatment to groups
    for ($GroupNo=0; $GroupNo<$TGCount; $GroupNo++) {
        //echo "\tAssigning values for treatment group $GroupNo\n";
        if (($GroupAttrs[$GroupNo]["Area"]>$MinAreaForAerial) &&
                ($GroupAttrs[$GroupNo]["Pct_Cov"]>$MinPctCovForAerial))
        {
            $GroupAttrs[$GroupNo]["Treatment"]="Aerial";
        } else {
            if (($GroupAttrs[$GroupNo]["SlopeMAX"]<=$MaxSlopeForMechanical) &&
                    ($GroupAttrs[$GroupNo]["Pct_Cov"]>=$MinPctCovForMechanical))
            {
                $GroupAttrs[$GroupNo]["Treatment"]="Mechanical";
            } else {
                $GroupAttrs[$GroupNo]["Treatment"]="Manual";
            }
        }
    }

    // Dis-aggregate treatment groups
    for ($TG=0; $TG<$TGCount; $TG++)
    {
        $NumPolysInTG=count($TreatmentGroups[$TG]);
        for ($TGPoly=0;$TGPoly<$NumPolysInTG;$TGPoly++)
        {
            $Poly=$TreatmentGroups[$TG][$TGPoly];
            $TamariskPolygons[$Poly]["Treatment"]=$GroupAttrs[$TG]["Treatment"];
        }

    }

    return($TamariskPolygons);
}

function ReadTamariskPolygons($TamariskCSVFile)
{
// Read the tamarisk data file and assign it to an array
    if (($handle = fopen("$TamariskCSVFile", "r")) == FALSE) {
        die("ReadTamariskPolygons: Failed to open file $TamariskCSVFile.");
    } else {
        echo "ReadTamariskPolygons: Opened file\n";
        // Get the column headers, look for columns we want to use
        $data = fgetcsv($handle, 10000, ",");
        $NumColumns = count($data);
        echo "ReadTamariskPolygons: $NumColumns columns in tamarisk file\n";

        for ($c=0; $c < $NumColumns; $c++) {
            $ColumnName=$data[$c];
            switch ($ColumnName) {
                case "VALUE":
                    $Value=$c;
                    break;
                case "Height":
                    $Height=$c;
                    break;
                case "Pct_Cov":
                    $Pct_Cov=$c;
                    break;
                case "AREA":
                    $Area=$c;
                    break;
                case "SlopeMAX":
                    $SlopeMAX=$c;
                    break;
                case "SlopeMEAN":
                    $SlopeMean=$c;
                    break;
                case "CostMEAN":
                    $CostMEAN=$c;
                   
                    break;
            }
        }
        echo "Columns:\n";
        echo "\tValue=$Value, Height=$Height, Pct_Cov=$Pct_Cov\n";
        echo "\tArea=$Area, SlopeMAX=$SlopeMAX, SlopeMean=$SlopeMean, CostMEAN=$CostMEAN\n";

        // Now get the data rows and store them in an array
        $TamariskPolygons=array();
        $row = 0;
        while (($data = fgetcsv($handle, 1000, ",")) !== FALSE) {
            for ($c=0; $c < $NumColumns; $c++) {
                //echo "Row: $row, ".$data[$c] . "<br />\n";
                $AttributeValue=$data[$c];
                switch ($c) {
                    case $Value:
                        $TamariskPolygons[$row]["Value"]=$AttributeValue;
                        break;
                    case $Height:
                        $TamariskPolygons[$row]["Height"]=$AttributeValue;
                        break;
                    case $Pct_Cov:
                        $TamariskPolygons[$row]["Pct_Cov"]=$AttributeValue;
                        break;
                    case $Area:
                        $TamariskPolygons[$row]["Area"]=$AttributeValue;
                        break;
                    case $SlopeMAX:
                        $TamariskPolygons[$row]["SlopeMAX"]=$AttributeValue;
                        break;
                    case $SlopeMean:
                        $TamariskPolygons[$row]["SlopeMean"]=$AttributeValue;
                        break;
                    case $CostMEAN:
                        $TamariskPolygons[$row]["CostMEAN"]=$AttributeValue;
                        break;
                }
            }
            // Calculate scaled cost, using a base of BASECOST units
            $TamariskPolygons[$row]["CostMultiplier"]=
                1+(($TamariskPolygons[$row]["CostMEAN"]-BASECOST)/BASECOST);
            $row++;
        }
        fclose($handle);
    }
    echo "There are $row polygons.\n";
    return($TamariskPolygons);
}

function ReadCosts($FileName)
// Costs are read in cost/acre. Convert to cost/square meter
// by dividing by 4046.85642
{
    if (($handle = fopen("$FileName", "r")) == FALSE) {
        die("ReadCosts: Failed to open file $FileName.");
    } else {
        echo "ReadCosts: Opened file $FileName\n";

        $Matrix=array();
        $TreatmentCount=0;

        // Discard the first row - it is the header row
        $data = fgetcsv($handle, 1000, ",");

        // Columns are: 0=treatment type, 1=treatment cost/acre
        while (($data = fgetcsv($handle, 1000, ",")) !== FALSE) {

            // Get the treatment name and the cost
            $Treatment=$data[0];
            $Cost=$data[1];
            $AdjustedCost=$Cost/4046.85642;
            $Matrix[$Treatment]=$AdjustedCost;
            echo "\tReadCosts: Assigned treatment $Treatment cost $AdjustedCost.\n";
            $TreatmentCount++;
        }
        fclose($handle);
        $NumTreatments=$TreatmentCount;
        echo "\tReadCosts: Found $NumTreatments treatments.\n";

        return($Matrix);
    }
}

function PrintGroup($GroupNo, $GroupMatrix)
{
    $NoPolysInGroup=count($GroupMatrix[$GroupNo]);
    //echo "$GroupNo: NoPolysInGroup=$NoPolysInGroup\n";
    for ($i=0;$i<$NoPolysInGroup;$i++)
    {
        $PolyNo=$GroupMatrix[$GroupNo][$i];
        echo "$PolyNo ";
    }
    echo "\n";
}

function AssignGroups ($TamariskPolygons, $DistanceMatrix, $OldGroups)
// returns a new trial grouping, based on distance between polygons,
// existing grouping, and preferred treatments
{
    // Aggregate treatment groups
    //echo "Calling AggregateTGs from AssignGroups\n";
    $GroupAttrs=AggregateTGs ($OldGroups, $TamariskPolygons);

    $NumPolygons=count($TamariskPolygons);
    $OldTGCount=count($OldGroups);
    $NewGroups=array();

    // We are going to merge polygons into groups, and merge groups.
    // Once a group or polygon is merged, it is not
    // available to be merged again. Set a flag  to indicate
    // if it has been merged.
    $GroupMergedFlag=array();
    $PolyMergedFlag=array();
    // Set a counter for each new group to indicate how many polygons are members
    // of that group. Since new groups have not yet been assigned, this number
    // starts off at zero. There can not be more new groups than old groups.
    $NewGroupPolyCount=array();
    // A couter for the number of polygons in the old groupings
    $OldGroupPolyCount=array();

    // Need to keep track of which groups have already been examined
    for ($i=0;$i<$OldTGCount;$i++)
    {
        $GroupMergedFlag[$i]=0;
        $NewGroupPolyCount[$i]=0;
        $OldGroupPolyCount[$i]=count($OldGroups[$i]);
    }
    // Also keep track of which polygons have been examined
    for ($i=0;$i<$NumPolygons;$i++)
    {
        $PolyMergedFlag[$i]=0;
    }

    // We want to consider Aerial treatment groups first, then mechanical,
    // and finally manual
    $GroupOrder[0]="Aerial";
    $GroupOrder[1]="Mechanical";
    $GroupOrder[2]="Manual";
    $NumTreatmentTypes=count($GroupOrder);

    // Count of the groups in the new grouping
    $NewGroupNo=0;

    // need to look at each group, and figure out if it can be joined
    // with the nearest group. The processing order will be from least expensive
    // treatment type to most expensive treatment type
    for ($TreatmentTypeNo=0;$TreatmentTypeNo<$NumTreatmentTypes;$TreatmentTypeNo++)
    {
        $TreatmentName=$GroupOrder[$TreatmentTypeNo];

        for ($OldGroupNo=0;$OldGroupNo<$OldTGCount;$OldGroupNo++)
        // If this group has not yet been merged with another group,
        // start a new group in the output grouping.
        // If this group has already been merged with another group,
        // do not consider it for further grouping
        {
            // Only process this group if it has the current treatment type
            if ((!$GroupMergedFlag[$OldGroupNo]) &&
                    ($TamariskPolygons[$OldGroups[$OldGroupNo][0]]["Treatment"]==$TreatmentName))
            {
                // Set the flag to show that this group has been processed
                $GroupMergedFlag[$OldGroupNo]=1;

                // Add the polygons currently in this group to the new grouping
                // and note that these polygons have been processed.
                for ($i=0; $i<$OldGroupPolyCount[$OldGroupNo]; $i++)
                {
                    $PolyNo=$OldGroups[$OldGroupNo][$i];
                    $NewGroups[$NewGroupNo][$i]=$PolyNo;
                    $PolyMergedFlag[$PolyNo]=1;
                }
                $NewGroupPolyCount[$NewGroupNo]=count($NewGroups[$NewGroupNo]);

                // Find the closest group to the current group
                $NearestGroup=CentroidMatrix::NearestGroup
                    ($OldGroupNo, $OldGroups, $TamariskPolygons, $DistanceMatrix);
                $SecondNearestGroup=CentroidMatrix::NthNeighborGroup
                    ($OldGroupNo, $OldGroups, $TamariskPolygons, $DistanceMatrix, 1);
                $NearestGroupDistance=CentroidMatrix::DistanceBetweenGroups
                    ($NearestGroup, $OldGroupNo, $DistanceMatrix, $OldGroups);
                $SecondNearestGroupDistance=CentroidMatrix::DistanceBetweenGroups
                    ($SecondNearestGroup, $OldGroupNo, $DistanceMatrix, $OldGroups);
/*
                echo "Considering group $OldGroupNo:";
                PrintGroup($OldGroupNo, $OldGroups);
                echo "\tNearest ($NearestGroup:$NearestGroupDistance): ";
                PrintGroup($NearestGroup, $OldGroups);
                echo "\tSecond Nearest ($SecondNearestGroup:$SecondNearestGroupDistance): ";
                PrintGroup($SecondNearestGroup, $OldGroups);
*/
                $GroupToMerge=-1;
                if ((!$GroupMergedFlag[$NearestGroup]) &&
                        (CanMergeGroups($OldGroupNo, $NearestGroup,
                                $GroupAttrs, $NearestGroupDistance)))
                {
                    $GroupToMerge=$NearestGroup;
                    echo "Merging $OldGroupNo to nearest group $GroupToMerge\n";
                } elseif ((!$GroupMergedFlag[$SecondNearestGroup]) &&
                        (CanMergeGroups($OldGroupNo, $SecondNearestGroup,
                                $GroupAttrs, $SecondNearestGroupDistance))) {
                    $GroupToMerge=$SecondNearestGroup;
                    echo "Merging $OldGroupNo to second nearest group $GroupToMerge\n";
                }
                if ($GroupToMerge>=0)
                {
                    //echo "We can merge groups $OldGroupNo and $GroupToMerge.\n";
                    // If we can merge the groups, Add the old group to the new group
                    for ($g=0;$g<$OldGroupPolyCount[$GroupToMerge]; $g++)
                    {
                        $NewGroupPolyNo=$NewGroupPolyCount[$NewGroupNo];
                        $PolyNo=$OldGroups[$GroupToMerge][$g];
                        $NewGroups[$NewGroupNo][$NewGroupPolyNo]=$PolyNo;
                        $PolyMergedFlag[$PolyNo]=1;
                        $NewGroupPolyCount[$NewGroupNo]++;
                    }

                    // This group has been merged; it is not available for
                    // further merges in this iteration
                    $GroupMergedFlag[$GroupToMerge]=1;
                } else {
                    // We didn't fing a group to merge. Try to "steal" a
                    // Polygon from an adjacent group.
                    $PolyDistanceArray=CentroidMatrix::RankPolysByDistance
                        ($OldGroupNo, $OldGroups, $TamariskPolygons,
                         $DistanceMatrix, 1);

                    // It would be better to define a radius to search from
                    // each polygon in this group, and continue searching until
                    // we are "too far" away to join polygons

                    // Try to steal the nearest poly to this group, or the second
                    // nearest if the fist is not OK
                    $OKToStealPoly=0;
                    $PolyToSteal=$PolyDistanceArray[0]["PolyNo"];
                    //echo "About to try to integrate(0) polygon $PolyToSteal into group $OldGroupNo\n";
                    if ((!$PolyMergedFlag[$PolyToSteal]) &&
                        CanMergePolyToGroup($OldGroupNo, $GroupAttrs,
                            $PolyToSteal, $TamariskPolygons, $DistanceMatrix))
                    {
                        $OKToStealPoly=1;
                    } else {
                        $PolyToSteal=$PolyDistanceArray[1]["PolyNo"];
                        //echo "About to try to integrate (1) polygon $PolyToSteal into group $OldGroupNo\n";
                        if ((!$PolyMergedFlag[$PolyToSteal]) &&
                            CanMergePolyToGroup($OldGroupNo, $GroupAttrs,
                                $PolyToSteal, $TamariskPolygons, $DistanceMatrix))
                        {
                            $OKToStealPoly=1;
                        }
                    }

                    if ($OKToStealPoly)
                    {
                        // Need to remove the polygon from the old group,
                        $ContainingGroup=GroupContainingPoly($OldGroups, $PolyToSteal);
                        echo "Moving polygon $PolyToSteal from group $ContainingGroup to group $OldGroupNo\n";
                        $TempGroup=array();
                        $TempPolyNo=0;
                        $NumPolysInGroup=count($OldGroups[$ContainingGroup]);
                        for ($P=0;$P<$NumPolysInGroup;$P++)
                        {
                            $PolyNo=$OldGroups[$ContainingGroup][$P];
                            if ($PolyNo != $PolyToSteal)
                            {
                                $TempGroup[$TempPolyNo]=$PolyNo;
                                $TempPolyNo++;
                            }
                        }
                        $OldGroups[$ContainingGroup]=$TempGroup;
                       
                        // It is possible to steal the only polygon from a group
                        // If so, remove the group from future consideration
                        if (count($OldGroups[$ContainingGroup]==0))
                            $GroupMergedFlag[$ContainingGroup]=1;

                        // Add polygon to the new group
                        $NewGroupPolyNo=$NewGroupPolyCount[$NewGroupNo];
                        $NewGroups[$NewGroupNo][$NewGroupPolyNo]=$PolyToSteal;
                        $NewGroupPolyCount[$NewGroupNo]++;

                        // Set the flag to show that this polygon has been merged
                        $PolyMergedFlag[$PolyToSteal]=1;
                    }

                }
                $NewGroupNo++;
            }
        }
    }
    return($NewGroups);
}

function CanMergeGroups($GroupNo1, $GroupNo2, $GroupAttrs, $Distance)
// Check if the attributes for these groups are similar enough to merge them
{
    global $MinPctCovForAerial;
    global $MinAreaForAerial;
    global $MaxSlopeForMechanical;
    global $MinPctCovForMechanical;

    $CanMergeGroups=1;

    if (($GroupAttrs[$GroupNo1]["Pct_Cov"]>=$MinPctCovForAerial) &&
            ($GroupAttrs[$GroupNo2]["Pct_Cov"]<$MinPctCovForAerial))
        $CanMergeGroups=0;
    if (($GroupAttrs[$GroupNo2]["Pct_Cov"]>=$MinPctCovForAerial) &&
            ($GroupAttrs[$GroupNo1]["Pct_Cov"]<$MinPctCovForAerial))
        $CanMergeGroups=0;

    if (($GroupAttrs[$GroupNo1]["Pct_Cov"]>=$MinPctCovForMechanical) &&
            ($GroupAttrs[$GroupNo2]["Pct_Cov"]<$MinPctCovForMechanical))
        $CanMergeGroups=0;
    if (($GroupAttrs[$GroupNo2]["Pct_Cov"]>=$MinPctCovForMechanical) &&
            ($GroupAttrs[$GroupNo1]["Pct_Cov"]<$MinPctCovForMechanical))
        $CanMergeGroups=0;

    if (($GroupAttrs[$GroupNo1]["SlopeMAX"]>=$MaxSlopeForMechanical) &&
            ($GroupAttrs[$GroupNo2]["SlopeMAX"]<$MaxSlopeForMechanical))
        $CanMergeGroups=0;
    if (($GroupAttrs[$GroupNo2]["SlopeMAX"]>=$MaxSlopeForMechanical) &&
            ($GroupAttrs[$GroupNo1]["SlopeMAX"]<$MaxSlopeForMechanical))
        $CanMergeGroups=0;

    //echo "CanMergeGroups: Groups $GroupNo1 and $GroupNo2 merge: $CanMergeGroups.\n";
    return($CanMergeGroups);
}

function CanMergePolyToGroup($GroupNo, $GroupAttrs, $PolyNo, $TamariskPolygons, $Distance)
// Check if the attributes for this polygon is similar enough to the group
// Currently $Distance is not used, but it may be used in the future
{
    global $MinPctCovForAerial;
    global $MinAreaForAerial;
    global $MaxSlopeForMechanical;
    global $MinPctCovForMechanical;

    $CanMergeGroups=1;

    if (($GroupAttrs[$GroupNo]["Pct_Cov"]>=$MinPctCovForAerial) &&
            ($TamariskPolygons[$PolyNo]["Pct_Cov"]<$MinPctCovForAerial))
        $CanMergeGroups=0;
    if (($TamariskPolygons[$PolyNo]["Pct_Cov"]>=$MinPctCovForAerial) &&
            ($GroupAttrs[$GroupNo]["Pct_Cov"]<$MinPctCovForAerial))
        $CanMergeGroups=0;

    if (($GroupAttrs[$GroupNo]["Pct_Cov"]>=$MinPctCovForMechanical) &&
            ($TamariskPolygons[$PolyNo]["Pct_Cov"]<$MinPctCovForMechanical))
        $CanMergeGroups=0;
    if (($TamariskPolygons[$PolyNo]["Pct_Cov"]>=$MinPctCovForMechanical) &&
            ($GroupAttrs[$GroupNo]["Pct_Cov"]<$MinPctCovForMechanical))
        $CanMergeGroups=0;

    if (($GroupAttrs[$GroupNo]["SlopeMAX"]>=$MaxSlopeForMechanical) &&
            ($TamariskPolygons[$PolyNo]["SlopeMAX"]<$MaxSlopeForMechanical))
        $CanMergeGroups=0;
    if (($TamariskPolygons[$PolyNo]["SlopeMAX"]>=$MaxSlopeForMechanical) &&
            ($GroupAttrs[$GroupNo]["SlopeMAX"]<$MaxSlopeForMechanical))
        $CanMergeGroups=0;

    //echo "CanMergePolyToGroup: Group $GroupNo and polygon $PolyNo merge: $CanMergeGroups.\n";
    return($CanMergeGroups);
}

function PrintTreatmentGroups ($TreatmentGroups, $TamariskPolygons)
{
    for ($TG=0;$TG<count($TreatmentGroups);$TG++)
    {
        echo "\tGroup $TG: ";
        for ($TGPoly=0; $TGPoly<count($TreatmentGroups[$TG]);$TGPoly++)
        {
            $Poly=$TreatmentGroups[$TG][$TGPoly];
            if ($TGPoly==0)
            {
                $Treatment=$TamariskPolygons[$Poly]["Treatment"];
                echo "$Treatment: ";
            }
            echo "$Poly ";
        }
        echo "\n";
    }
}

function GroupContainingPoly($Groups, $PolyToFind)
{
    //echo "GroupContainingPoly: Looking for $PolyToFind\n";
    $ContainingGroup=-1;
    $GroupCount=count($Groups);
    //echo "GroupContainingPoly: Number of groups to search is $GroupCount\n";
    for ($TG=0; $TG<$GroupCount; $TG++)
    {
        $PolyCount=count($Groups[$TG]);
        //echo "GroupContainingPoly: Group $TG has $PolyCount polys\n";
        for ($Poly=0;$Poly<$PolyCount;$Poly++)
        {
            $PolyNo=$Groups[$TG][$Poly];
            if ($PolyNo==$PolyToFind)
                $ContainingGroup=$TG;
        }
    }
    return($ContainingGroup);
}
////////////////////////////////////////////////////////////////////////////
//End of Functions
////////////////////////////////////////////////////////////////////////////

// Input files
//$TamariskCSVFile="data/NAD1983HARNColoradoSouth0503/WildlifeRefuge/Tamarisk.csv";
//$CostCSVFile="data/NAD1983HARNColoradoSouth0503/WildlifeRefuge/TreatmentCosts.csv";
//$CentroidMatrixCSVFile="data/NAD1983HARNColoradoSouth0503/WildlifeRefuge/CentroidMatrix.csv";
//$OutputFile = "data/NAD1983HARNColoradoSouth0503/WildlifeRefuge/Prescription.csv";
$TamariskCSVFile="data/NAD1983HARNColoradoSouth0503/WesternCO/Tamarisk.csv";
$CostCSVFile="data/NAD1983HARNColoradoSouth0503/WesternCO/TreatmentCosts.csv";
$CentroidMatrixCSVFile="data/NAD1983HARNColoradoSouth0503/WesternCO/CentroidMatrix.csv";
$OutputFile = "data/NAD1983HARNColoradoSouth0503/WesternCO/Prescription.csv";

///////////////////////////////////////
// Read in Tamarisk file
///////////////////////////////////////
$TamariskPolygons=ReadTamariskPolygons($TamariskCSVFile);
$NumPolygons=count($TamariskPolygons);
//echo "Number of polygons: $NumPolygons\n";
//print_r($TamariskPolygons);

///////////////////////////////////////
// Read in centroid matrix file
///////////////////////////////////////
$DistanceMatrix=CentroidMatrix::ReadCSVFile($CentroidMatrixCSVFile);
//print_r($DistanceMatrix);

///////////////////////////////////////
// Read in treatment cost file
///////////////////////////////////////
$TreatmentCosts=ReadCosts($CostCSVFile);
//print_r($TreatmentCosts);

///////////////////////////////////////
// Assign treatment groups
// Initially every polygon is in a separate treatment group
///////////////////////////////////////
$ProposedTreatmentGroups=array();
$TreatmentGroups=array();
for ($i=0;$i<$NumPolygons;$i++)
{
    $ProposedTreatmentGroups[$i][0]=$i;
    $TamariskPolygons[$i]["Group"]=$i;
}

$NewScaledTotalCost=0;
$NewTotalCost=0;
$ScaledTotalCost=1;
$TotalCost=1;
$Trials=1;
while (($NewScaledTotalCost<=$ScaledTotalCost) && ($Trials<200)
        && ($ProposedTreatmentGroups!= $TreatmentGroups))
{
    $ScaledTotalCost=$NewScaledTotalCost;
    $TotalCost=$NewTotalCost;

    //////////////////////////////////////
    // Assign a treatment
    //////////////////////////////////////
    $NewTamariskPolygons=AssignTreatment($TamariskPolygons,
            $ProposedTreatmentGroups);
    echo "Here are the suggested Treatments:\n";
    PrintTreatmentGroups ($ProposedTreatmentGroups, $NewTamariskPolygons);

    // Calculate total treatment cost
    $NewTotalCost=0;
    $NewScaledTotalCost=0;
    for ($Poly=0;$Poly<$NumPolygons;$Poly++)
    {
        $Treatment=$NewTamariskPolygons[$Poly]["Treatment"];
        $Area=$NewTamariskPolygons[$Poly]["Area"];
        $TreatmentCost=$TreatmentCosts[$Treatment];
        $CostMultiplier=$NewTamariskPolygons[$Poly]["CostMultiplier"];
        $NewTotalCost+=$TreatmentCost*$Area;
        $NewScaledTotalCost+=$TreatmentCost*$Area*$CostMultiplier;
    }
    echo "Unajusted treatment cost is $NewTotalCost.\n";
    echo "Scaled treatment cost is $NewScaledTotalCost.\n";

    // Special case for first iteration
    if ($ScaledTotalCost==0) $ScaledTotalCost=$NewScaledTotalCost+1;

    // If the new grouping is better than the old grouping, try yet
    // another grouping
    if ($NewScaledTotalCost<=$ScaledTotalCost)
    {
        $TreatmentGroups=$ProposedTreatmentGroups;
        $TamariskPolygons=$NewTamariskPolygons;

        // Put the tamarisk polygons in the new treatment groups
        for ($TG=0; $TG<count($TreatmentGroups);$TG++)
        {
            for ($Poly=0;$Poly<count($TreatmentGroups[$TG]);$Poly++)
            {
                $TamariskPolygons[$TreatmentGroups[$TG][$Poly]]["Group"]=$TG;
            }
        }
        //print_r($TamariskPolygons);

        $ProposedTreatmentGroups=
            AssignGroups($TamariskPolygons, $DistanceMatrix, $TreatmentGroups);
        //echo "Proposed treatment groups:\n";
        //print_r($ProposedTreatmentGroups);
    }
    $Trials++;
}

echo "The final proposed treatment groups:\n";
PrintTreatmentGroups ($TreatmentGroups, $TamariskPolygons);
echo "Unajusted treatment cost is $TotalCost.\n";
echo "Scaled treatment cost is $ScaledTotalCost.\n";
echo "Used $Trials trials to find this solution.\n";

// Write results to a csv file
$fh = fopen($OutputFile, 'w') or die("can't open file $myFile");
fwrite($fh, "VALUE,Treatment,Group\n");
for ($Polygon=0; $Polygon<$NumPolygons; $Polygon++) {
    $PolyNumber=$TamariskPolygons[$Polygon]["Value"];
    $Treatment =$TamariskPolygons[$Polygon]["Treatment"];
    $Group =$TamariskPolygons[$Polygon]["Group"];
    $StringData="$PolyNumber,$Treatment,$Group\n";
    fwrite($fh, $StringData);
}
fclose($fh);

?>




The source code for CentroidMatrix.php:
<?php
// Class to manipulate a matrix of distances between polygons
//
// Ted Manahan
// NR 505 Fall 2010
//error_reporting(E_ALL);
ini_set('display_errors', '1');

// Define a maximum distance to use as an initialization value
define("MAXDIST",999999999);

class CentroidMatrix {

    // function to read the .csv file and create the matrix
    function ReadCSVFile ($FileName)
    {
        if (($handle = fopen("$FileName", "r")) == FALSE) {
            die("Failed to open file $FileName.");
        } else {
            $PrintData=0;
            echo "ReadCSVFile: Opened centroid file $FileName\n";

            // How big is the matrix
            $data = fgetcsv($handle, 1000, ",");
            $num = count($data);
            echo "ReadCSVFile: $num columns in centroid file\n";

            // look for missing columns. For some reason ArcGIS will fail to
            // generate the centroid for some polygons
            $ColumnCount=0;
            $HaveData=array();
            for ($i=1;$i<$num;$i++)
            {
                if ($PrintData) echo "Line $i: ";
                $CentroidNo=$data[$i];
                while ($CentroidNo>$ColumnCount)
                {
                    // Missing a centroid data point
                    if ($PrintData) echo "Missing $ColumnCount ";
                    $HaveData[$ColumnCount]=0;
                    $ColumnCount++;
                }
                $HaveData[$ColumnCount]=1;
                $ColumnCount++;
                if ($PrintData) echo "$CentroidNo, \n";
            }
            if ($PrintData) print_r($HaveData);

            $Matrix=array();

            // Data starts in the second row
            $RowNumber=0;
            while (($data = fgetcsv($handle, 10000, ",")) !== FALSE) {
                //print_r($data);
                while (!$HaveData[$RowNumber])
                {
                    if ($PrintData) echo "No data for row $RowNumber: ";
                    for ($ColumnNumber=0; $ColumnNumber < $ColumnCount; $ColumnNumber++)
                    {
                        if ($RowNumber==$ColumnNumber)
                        {
                            $Matrix[$RowNumber][$ColumnNumber]=0;
                            if ($PrintData) echo "0,";
                        } else {
                            $Matrix[$RowNumber][$ColumnNumber]=MAXDIST;
                            if ($PrintData) echo MAXDIST.",";
                        }
                    }
                    $RowNumber++;
                    if ($PrintData) echo "\n";
                }

                if ($PrintData) echo "Have data for row $RowNumber: ";
                $ColumnNumber=0;

                // Get the distance to each object from the current object
                for ($col=1; $col < $num; $col++) {
                    // Fill in missing values
                    while (!$HaveData[$ColumnNumber])
                    {
                        $Matrix[$RowNumber][$ColumnNumber]=MAXDIST;
                        if ($PrintData) echo "$ColumnNumber:999999999,";
                        $ColumnNumber++;
                    }
                    $AttributeValue=$data[$col];
                    if ($RowNumber==$ColumnNumber)
                    {
                        $Matrix[$RowNumber][$ColumnNumber]=0;
                        if ($PrintData) echo "$ColumnNumber:0,";
                    } else {
                        $Matrix[$RowNumber][$ColumnNumber]=$AttributeValue;
                        if ($PrintData) echo "$ColumnNumber:$AttributeValue,";
                    }
                    $ColumnNumber++;
                }
                if ($PrintData) echo "\n";
                $RowNumber++;
            }
            fclose($handle);
            //print_r($Matrix);
            return($Matrix);
        }
    }

    function NearestGroup($GroupNumber, $Groups, $TamariskPolygons,
            $DistanceMatrix)
    // Return the group that is closest to GroupNumber
    {
        $TGCount=count($Groups);
        $PolyCount=count($TamariskPolygons);

        //echo "NearestGroup: Looking for the closest group to group # $GroupNumber\n";
        $NumPolysInGroup=count($Groups[$GroupNumber]);
        $ClosestPolyDist=MAXDIST+1;
        $CloseList=array();

        // For each polygon in this treatment group
        for ($j=0; $j<$NumPolysInGroup; $j++)
        {
            $CloseList[$j]=-1;
            $PolyNo=$Groups[$GroupNumber][$j];

            // Look for the closest polygon
            for ($k=0; $k<$PolyCount; $k++)
            {
                $Distance=$DistanceMatrix[$PolyNo][$k];
                if (($TamariskPolygons[$k]["Group"]<>$GroupNumber) &&
                    ($Distance<$ClosestPolyDist) &&
                    ($TamariskPolygons[$PolyNo]["Value"]<>$TamariskPolygons[$k]["Value"]))
                {
                    $ClosestPolyDist=$Distance;
                    $ClosestPoly=$k;
                }
            }
        }
        $ClosestGroup=$TamariskPolygons[$ClosestPoly]["Group"];
        //echo "NearestGroup: Found $ClosestGroup\n";
        return($ClosestGroup);
    }

    function NthNeighborGroup($GroupNumber, $Groups, $TamariskPolygons,
            $DistanceMatrix, $N=0)
    // Return the group that is the Nth closest to GroupNumber
    {
        $TGCount=count($Groups);
        $PolyCount=count($TamariskPolygons);

        //echo "NthNeighbor: Looking for the $N neighbor to group # $GroupNumber\n";
        $NumPolysInTargetGroup=count($Groups[$GroupNumber]);
        $CloseList=array();

        // Counter of the treatment group we are considering
        $GroupCount=0;

        // For each treatment group
        for ($TG=0; $TG<$TGCount; $TG++)
        {
            if ($TG<>$GroupNumber)
            {
                //echo "Looking at TG=$TG\n";
                $CloseList[$GroupCount]["GroupNo"]=-1;
                $CloseList[$GroupCount]["Distance"]=MAXDIST+1;
                $NumPolysInCantidateGroup=count($Groups[$TG]);

                // Look at each polygon in the cantidate treatment group
                for ($k=0; $k<$NumPolysInCantidateGroup; $k++)
                {
                    $CantidatePolyNo=$Groups[$TG][$k];
                    //echo "CantidatePolyNo=$CantidatePolyNo\n";
                    $ClosestPolyDist=MAXDIST+1;
                    // and see how far it is from each polygon in the target group
                    for ($P=0;$P<$NumPolysInTargetGroup; $P++)
                    {
                        $TargetPolyNo=$Groups[$GroupNumber][$P];
                        //echo "TargetPolyNo=$TargetPolyNo\n";
                        $Distance=$DistanceMatrix[$CantidatePolyNo][$TargetPolyNo];
                        //echo "Distance matrix extry $CantidatePolyNo, $TargetPolyNo is $Distance\n";
                        if ($Distance<$ClosestPolyDist)
                        {
                            $ClosestPolyDist=$Distance;
                            $ClosestPoly=$TargetPolyNo;
                        }
                    }
                }

                // Add the new polygon distance and group info to the list
                $j=$GroupCount;
                while ($j>0 && $CloseList[$j-1]["Distance"]>$ClosestPolyDist)
                {
                    $CloseList[$j]["GroupNo"]=$CloseList[$j-1]["GroupNo"];
                    $CloseList[$j]["Distance"]=$CloseList[$j-1]["Distance"];
                    $j--;
                }
                $CloseList[$j]["GroupNo"]=$TG;
                $CloseList[$j]["Distance"]=$ClosestPolyDist;

                $GroupCount++;
            }
        }
        //print_r($CloseList);
        $N_CloseGroup=$CloseList[$N]["GroupNo"];
        $N_CloseDist =$CloseList[$N]["Distance"];
        //echo "NthNeighborGroup: $N neighbor to group $GroupNumber is ".
        //    "$N_CloseGroup with distance $N_CloseDist\n";
        return($CloseList[$N]["GroupNo"]);
     }

    function PolyInGroup($PolyNo, $GroupNo, $Groups)
    // Return 1 if the polygon is in this group, else 0
    {
        $IsInGroup=0;
        $NumPolysInGroup=count($Groups[$GroupNo]);

        for ($i=0;$i<$NumPolysInGroup;$i++)
        {
            if ($Groups[$GroupNo][$i]==$PolyNo)
                $IsInGroup=1;
        }
        //echo "PolyInGroup: Polygon $PolyNo in group $GroupNo: $IsInGroup\n";
        return($IsInGroup);
    }

    function RankPolysByDistance($GroupNumber, $Groups, $TamariskPolygons,
            $DistanceMatrix, $N=0)
    // Return an array of Polygons ranked by; distance from $Groups[$GroupNumber]
    {
        $TGCount=count($Groups);
        $PolyCount=count($TamariskPolygons);
        $PolysInCurrentGroup=count($Groups[$GroupNumber]);

        //echo "NthNeighborPoly: Looking for the $N neighbor to group # $GroupNumber\n";

        $CloseList=array();

        // Counter for the polygon we are considering
        $PolyNumber=0;
        $CloseList[$PolyNumber]["PolyNo"]=-1;
        $CloseList[$PolyNumber]["Distance"]=MAXDIST+1;

        // For each polygon
        for ($Poly=0; $Poly<$PolyCount; $Poly++)
        {
            if (!CentroidMatrix::PolyInGroup($Poly, $GroupNumber, $Groups))
            {
                //echo "NthNeighborPoly: Looking at Poly=$Poly\n";
                $ClosestPolyDist=MAXDIST+1;

                // See how far it is from the polygon to our group
                for ($P=0;$P<$PolysInCurrentGroup; $P++)
                {
                    $TargetPolyNo=$Groups[$GroupNumber][$P];
                    $Distance=$DistanceMatrix[$TargetPolyNo][$Poly];
                    //echo "\tNthNeighborPoly: Distance matrix extry $Poly, $TargetPolyNo is $Distance\n";
                    if ($Distance<$ClosestPolyDist)
                    {
                        $ClosestPolyDist=$Distance;
                        $ClosestPoly=$Poly;
                    }
                }

                // Add the new polygon distance and group info to the list
                $j=$PolyNumber;
                while ($j>0 && $CloseList[$j-1]["Distance"]>$ClosestPolyDist)
                {
                    $CloseList[$j]["PolyNo"]=$CloseList[$j-1]["PolyNo"];
                    $CloseList[$j]["Distance"]=$CloseList[$j-1]["Distance"];
                    $j--;
                }
                $CloseList[$j]["PolyNo"]=$ClosestPoly;
                $CloseList[$j]["Distance"]=$ClosestPolyDist;

                $PolyNumber++;
            }
        }
        //print_r($CloseList);
        return($CloseList);
     }

    function DistanceBetweenGroups ($Group1, $Group2, $DistanceMatrix, $Groups)
    {
        $G1Count=count($Groups[$Group1]);
        $G2Count=count($Groups[$Group2]);
        $Distance=MAXDIST+1;
        for ($i=0;$i<$G1Count;$i++)
        {
            $P1=$Groups[$Group1][$i];
            for ($j=0;$j<$G2Count;$j++)
            {
                $P2=$Groups[$Group2][$j];
                if ($DistanceMatrix[$P1][$P2]<$Distance)
                    $Distance=$DistanceMatrix[$P1][$P2];
            }
        }
        return($Distance);
    }
}
?>