Landscape.cpp

Go to the documentation of this file.

/*
*******************************************************************************************************
Copyright (c) 2011, Christopher John Topping, University of Aarhus
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided
that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the
following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
********************************************************************************************************
*/
#include <iostream>
#include <fstream>
#include <string>
#include "../Landscape/ls.h"
#include "../BatchALMaSS/PopulationManager.h"
#include "../Skylark/skylarks_all.h"
#include "../RodenticideModelling/RodenticidePredators.h"
#include "../Landscape/map_cfg.h"
#include "../BatchALMaSS/BoostRandomGenerators.h"


using namespace std;
/*
This is where all the boost random number generators are created
*/
base_generator_type g_generator(static_cast<unsigned int>(std::time(0)));
boost::uniform_real<> g_uni_dist(0,1);
boost::uniform_int<> g_uni_dist2(0,10000);
boost::variate_generator<base_generator_type&, boost::uniform_real<> > g_rand_uni(g_generator, g_uni_dist);
boost::variate_generator<base_generator_type&, boost::uniform_int<> > g_rand_uni2(g_generator, g_uni_dist2);
// End boost random number generators

// Version info.
static const int version_major = 1;
static const int version_minor = 0;
static const int version_revision = 0;
static const char * version_date = "2014-02-03";


extern CfgBool l_pest_enable_pesticide_engine;
CfgBool cfg_rodenticide_enable("RODENTICIDE_ENABLE",CFG_CUSTOM,false);
CfgBool cfg_rodenticide_reporting_enable("RODENTICIDE_REPORTING_ENABLE",CFG_CUSTOM,false);

static CfgBool cfg_map_usesoiltypes("MAP_USESOILTYPES",CFG_CUSTOM,false);
static CfgBool cfg_CalculateCentroids("MAP_CALCULATE_CENTROIDS",CFG_CUSTOM,false);

extern CfgBool cfg_OptimiseBedriftsmodelCrops;
extern CfgBool cfg_OptimisingFarms;
extern CfgBool cfg_DumpFarmAreas;
extern CfgBool cfg_MaizeEnergy;

// Hedge Subtype controls
CfgInt cfg_HedgeSubtypeMinimum("HEDGE_SUBTYPEMINIMUM",CFG_CUSTOM,0);
CfgInt cfg_HedgeSubtypeMaximum("HEDGE_SUBTYPEMAXIMUM",CFG_CUSTOM,3); // actual max is 2 (<3)

// Beetlebank variables
static CfgBool cfg_AddBeetleBanks("BBANKS_ADD",CFG_CUSTOM,false);
static CfgInt cfg_BeetleBankWidth("BBANK_WIDTH",CFG_CUSTOM,4);
static CfgInt cfg_BeetleBankChance("BBANK_CHANCE",CFG_CUSTOM,100); // out of 100
static CfgFloat cfg_BeetleBankMaxArea("BBANK_MAXAREA",CFG_CUSTOM,0.05);
CfgBool cfg_BeetleBankInvert("BBANK_INVERT",CFG_CUSTOM,false);
CfgInt cfg_BeetleBankMinX("BBANK_MINX",CFG_CUSTOM,0);
CfgInt cfg_BeetleBankMinY("BBANK_MINY",CFG_CUSTOM,0);
CfgInt cfg_BeetleBankMaxX("BBANK_MAXX",CFG_CUSTOM,100000);
CfgInt cfg_BeetleBankMaxY("BBANK_MAXY",CFG_CUSTOM,100000);
static CfgInt cfg_BeetleBankType("BBANK_TYPE",CFG_CUSTOM, (int)tole_BeetleBank);

// Veg area dump config variables
CfgBool cfg_dumpvegjan( "G_VEGAREASJAN_ON", CFG_CUSTOM, false );
CfgStr cfg_dumpvegjanfile( "G_VEGAREASJAN_FILENAME", CFG_CUSTOM, "DUMPVEG_JAN.TXT" );
CfgBool cfg_dumpvegjune( "G_VEGAREASJUNE_ON", CFG_CUSTOM, false );
CfgStr cfg_dumpvegjunefile( "G_VEGAREASJUNE_FILENAME", CFG_CUSTOM, "DUMPVEG_JUNE.TXT" );

// Crops configuration options
CfgFloat cfg_strigling_prop( "CROPS_STRIGLING_PROPORTION", CFG_CUSTOM, 1.0 );
CfgFloat cfg_silage_prop( "CROPS_SILAGE_PROPORTION", CFG_CUSTOM, 1.0 );
CfgFloat cfg_ins_app_prop1("CROPS_INSECTICIDE_APPLIC_ONE_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_ins_app_prop2("CROPS_INSECTICIDE_APPLIC_TWO_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_ins_app_prop3("CROPS_INSECTICIDE_APPLIC_THREE_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_herbi_app_prop("CROPS_HERBICIDE_APPLIC_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_fungi_app_prop1("CROPS_FUNGICIDE_APPLIC_ONE_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_fungi_app_prop2("CROPS_FUNGICIDE_APPLIC_TWO_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_fungi_app_prop3("CROPS_FUNGICIDE_APPLIC_THREE_PROPORTION", CFG_CUSTOM, 1.0);
CfgFloat cfg_greg_app_prop("CROPS_GROWTHREGULATOR_APPLIC_PROPORTION", CFG_CUSTOM, 1.0);

// Below define the area for custom actions, if valid_x, and valid_y fall in the
// area defined by x,y,size, then something can be set on that polygon
// The something is hard-coded each time it is used
CfgInt cfg_l_treatment_x( "LAND_TREATMENTX", CFG_CUSTOM, 0 );
CfgInt cfg_l_treatment_y( "LAND_TREATMENTY", CFG_CUSTOM, 0 );
CfgInt cfg_l_treatment_size( "LAND_TREATMENTSIZE", CFG_CUSTOM, -1 );
CfgBool cfg_l_usecustompoly( "LAND_USECUSTOMPOLY", CFG_CUSTOM, false );

// The config variable below notes the year from the start of the sim when
// the pesticide engine will start to spray its pesticide
// A test for the value m_toxShouldSpray in landscape is required before sprays
// can be used.
CfgInt cfg_productapplicstartyear("PEST_PROCTAPPLICSTARTYEAR",CFG_CUSTOM,9999999);
CfgInt cfg_productapplicendyear("PEST_PROCTAPPLICENDYEAR",CFG_CUSTOM,-1);
CfgInt cfg_pesticidetesttype("PESTICIDETESTYPE",CFG_CUSTOM, -1);

// Local configuration options:
static CfgInt cfg_MaxPondSize( "MAP_MAXPONDSIZE", CFG_CUSTOM, 5000 );
static CfgInt l_map_no_pesticide_fields( "MAP_NO_PESTICIDE_FIELDS", CFG_CUSTOM, 0 );
static CfgBool l_map_print_version_info( "MAP_PRINT_VERSION_INFO", CFG_CUSTOM, true );

// Whether to exit on the condition that the number of covered squares
// on the map becomes zero upon adding the border to a very slim polygon.
static CfgBool l_map_exit_on_zero_area( "MAP_EXIT_ON_ZERO_AREA", CFG_CUSTOM, true );

// Whether to check in the map for the existence of polygons that are
// mentioned in the polygon reference file.
static CfgBool l_map_check_polygon_xref( "MAP_CHECK_POLYGON_XREF", CFG_CUSTOM, true );

// Whether to add artificial hedgebanks to every hedge polygon in the map.
static CfgBool l_map_art_hedgebanks( "MAP_ART_HEDGEBANKS", CFG_CUSTOM, false );

static CfgStr l_map_map_file( "MAP_MAP_FILE", CFG_CUSTOM, "map.lsb" );
static CfgStr l_map_poly_file( "MAP_POLY_FILE", CFG_CUSTOM, "polygonrefs.txt" );
static CfgStr l_map_weather_file( "MAP_WEATHER_FILE", CFG_CUSTOM, "weather.pre" );
static CfgStr l_map_cropcurves_file( "MAP_CROPCURVES_FILE", CFG_CUSTOM, "curves.pre" );

static CfgInt l_map_chameleon_replace_num( "MAP_CHAMELEON_REPLACE_NUM", CFG_CUSTOM, 58 );

//'''''''''''''''''' CIPE LANDSCAPE MAKER CODE HERE  //'''''''''''''''''''''
static CfgBool l_map_CIPEmaker_enable( "MAP_CIPEMAKER_ENABLE", CFG_CUSTOM, false );
//'''''''''''''''''' CIPE LANDSCAPE MAKER CODE ABOVE  //'''''''''''''''''''''

static CfgBool l_map_dump_enable( "MAP_DUMP_ENABLE", CFG_CUSTOM, false );
static CfgBool l_map_removesmallpolygons("MAP_REMOVESMALLPOLYGONS", CFG_CUSTOM, false);
CfgStr l_map_dump_map_file( "MAP_DUMP_MAP_FILE", CFG_CUSTOM, "dump.lsb" );
static CfgBool l_map_dump_gfx_enable( "MAP_DUMP_GFX_ENABLE", CFG_CUSTOM, false );
static CfgStr l_map_dump_gfx_file( "MAP_DUMP_GFX_FILE", CFG_CUSTOM, "dump.ppm" );
static CfgBool l_map_dump_exit( "MAP_DUMP_EXIT", CFG_CUSTOM, false );
CfgStr l_map_dump_poly_file( "MAP_DUMP_POLY_FILE", CFG_CUSTOM, "dumppolyrefs.txt" );
static CfgInt l_map_umargin_width( "MAP_UMARGINWIDTH", CFG_CUSTOM, 12 );
static CfgStr l_map_dump_margin_file( "MAP_DUMP_MARGIN_FILE", CFG_CUSTOM, "dumpunsprayedmargins.txt" );
static CfgBool l_map_dump_treatcounts_enable( "MAP_DUMP_TREATCOUNTS_ENABLE", CFG_CUSTOM, false );
static CfgStr l_map_dump_treatcounts_file( "MAP_DUMP_TREATCOUNTS_FILE", CFG_CUSTOM, "treatment_counts.txt" );

static CfgBool l_map_dump_veg_enable( "MAP_DUMP_VEG_ENABLE", CFG_CUSTOM, true );
static CfgInt l_map_dump_veg_x( "MAP_DUMP_VEG_X", CFG_CUSTOM, 100 );
static CfgInt l_map_dump_veg_y( "MAP_DUMP_VEG_Y", CFG_CUSTOM, 100 );
static CfgBool l_map_dump_event_enable( "MAP_DUMP_EVENT_ENABLE", CFG_CUSTOM, false );
static CfgInt l_map_dump_event_x1( "MAP_DUMP_EVENT_XA", CFG_CUSTOM, 4287 );
static CfgInt l_map_dump_event_y1( "MAP_DUMP_EVENT_YA", CFG_CUSTOM, 2909 );
static CfgInt l_map_dump_event_x2( "MAP_DUMP_EVENT_XB", CFG_CUSTOM, 4333 );
static CfgInt l_map_dump_event_y2( "MAP_DUMP_EVENT_YB", CFG_CUSTOM, 2889 );
/*
static CfgBool l_map_renumberpolys("MAP_RENUMBERPOLY", CFG_CUSTOM, false);
*/
static CfgBool l_map_consolidatepolys("MAP_CONSOLIDATEPOLYS", CFG_CUSTOM, false);
static CfgBool l_map_calc_openness("MAP_CALC_OPENNESS",CFG_CUSTOM,false );
static CfgBool l_map_write_ascii( "MAP_WRITE_ASCII", CFG_CUSTOM, false );
static CfgInt l_map_ascii_utm_x( "MAP_ASCII_UTM_X", CFG_CUSTOM, 0 );
static CfgInt l_map_ascii_utm_y( "MAP_ASCII_UTM_Y", CFG_CUSTOM, 0 );
static CfgInt l_map_ascii_map_entity( "MAP_ASCII_MAP_ENTITY", CFG_CUSTOM, 1 );
/*
static CfgBool l_map_read_openness("MAP_READ_OPENNESS",CFG_CUSTOM,true );
static CfgBool l_map_write_openness("MAP_WRITE_OPENNESS",CFG_CUSTOM,false );
static CfgStr l_map_opennessfile("MAP_OPENNESSFILE",CFG_CUSTOM,"FieldOpennessScores.txt");
*/
// For building centroid calculations
static CfgInt cfg_mintownbuildingdistance("MAP_MINTOWNBUILDINGDISTANCE",CFG_CUSTOM, 100);
static CfgInt cfg_mintownbuildingnumber("MAP_MINTOWNBUILDINGNUMBER",CFG_CUSTOM, 6);

extern CfgBool cfg_rectangularmaps_on;

extern CfgFloat cfg_P1A;// Parameter 1
extern CfgFloat cfg_P1B; // Parameter 2
extern CfgFloat cfg_P1C; // Parameter 3
extern CfgFloat cfg_P1D; // The scaler to change to from g/dw to KJ
extern CfgBool  cfg_P1E; // reverse axis values
extern CfgFloat cfg_P1F; // Max x-value - at this point the curve tends to 0, must stope here to avoid negative values
extern CfgFloat cfg_P1G; // Min x-value
extern CfgStr   cfg_P1H;

// For hareas
extern double g_FarmIntensivenessH;
extern double g_VegHeightForageReduction;

extern CfgFloat cfg_G6A;// Parameter 1
extern CfgFloat cfg_G6B; // Parameter 2
extern CfgFloat cfg_G6C; // Parameter 3
extern CfgFloat cfg_G6D; // The scaler to change to from g/dw to KJ
extern CfgBool  cfg_G6E; // reverse axis values
extern CfgFloat cfg_G6F; // Max x-value - at this point the curve tends to 0, must stope here to avoid negative values
extern CfgFloat cfg_G6G; // Min x-value
extern CfgStr   cfg_G6H;

extern CfgFloat cfg_B6A;// Parameter 1
extern CfgFloat cfg_B6B; // Parameter 2
extern CfgFloat cfg_B6C; // Parameter 3
extern CfgFloat cfg_B6D; // The scaler to change to from g/dw to KJ
extern CfgBool  cfg_B6E; // reverse axis values
extern CfgFloat cfg_B6F; // Max x-value - at this point the curve tends to 0, must stope here to avoid negative values
extern CfgFloat cfg_B6G; // Min x-value
extern CfgStr   cfg_B6H;

extern CfgBool cfg_WriteCurve;



// Globals ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

double g_SpeedyDivides[2001];

int * m_polymapping;
Landscape * g_landscape_p;

// End Globals ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*
void Landscape::ReadOpenness()
{
/*
    int sz = (int) m_elems.size();
    ifstream ifile(l_map_opennessfile.value(),ios::in);
    if ( !ifile.is_open() )
    {
        Warn("Landscape::ReadOpenness() Cannot open input file ", "FieldOpennessScores.txt" );
    }
    int n, openness, polyref, cx, cy;
    ifile >> n;
    if (n != sz)
    {
        Warn("Landscape::ReadOpenness() number of entries differs from number of polygons ", l_map_opennessfile.value());
        std::exit(0);
    }
    for (int i=0; i<n; i++)
    {
        ifile >> polyref >> openness >> cx >> cy;
        // We have read a polygon reference in, but we need to find the associated m_elems entry. 
        // This can be done using the polymapping.
        int polyindex = m_polymapping[polyref];
        if (polyref != m_elems[polyindex]->GetPoly())
        {
            g_msg->Warn("Landscape::ReadOpenness() mismatch in polyref number at entry ", i);
            std::exit(0);
        }
        // All OK so save this
        m_elems[polyindex]->SetOpenness(openness);
        m_elems[polyindex]->SetCentroid(cx, cy);
    }
    ifile.close();
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::WriteOpenness()
{
/*
    int sz = (int) m_elems.size();
    ofstream ofile(l_map_opennessfile.value(),ios::out);
    if ( !ofile.is_open() )
    {
        Warn("Landscape::WriteOpenness() Cannot open output file ", l_map_opennessfile.value() );
    }
    ofile << sz << endl;
    for (int i=0; i<sz; i++)
    {
        ofile << m_elems[i]->GetPoly() << '\t' << m_elems[i]->GetOpenness() << '\t' << m_elems[i]->GetCentroidX()<<'\t'<<m_elems[i]->GetCentroidY()<<'\t'<< endl;
    }
    ofile.close();
}
//-----------------------------------------------------------------------------------------------------------------------------
*/
void Landscape::CalculateOpenness( bool a_realcalc)
{
    for (unsigned int i = 0; i < m_elems.size(); i++)
    {
        TTypesOfLandscapeElement tole = m_elems[i]->GetElementType();
        switch (tole)
        {
        case tole_Field:
        case tole_Marsh:
        case tole_Scrub:
        case tole_PermPastureLowYield:
        case tole_PermPastureTussocky:
        case tole_PermanentSetaside:
        case tole_PermPasture:
        case tole_NaturalGrassDry:
        case tole_NaturalGrassWet:
            if (a_realcalc)
            {
                cout << i << " ";
                m_elems[i]->SetOpenness(CalulateFieldOpennessAllCells(i));
            }
            else m_elems[i]->SetOpenness(0);
            break;
        default:
            m_elems[i]->SetOpenness(0);
            break;
        }
    }
    if (a_realcalc) cout << endl;
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::CalulateFieldOpennessCentroid(int  a_pref)
{
    // Get the centre point for this field
    int d0 = m_maxextent; // this is the width of the landscape and will always be at least as big as the biggest return value possible
    int cx = m_elems[a_pref]->GetCentroidX();
    int cy = m_elems[a_pref]->GetCentroidY();

    double offsetx = -1;
    double offsety = -1;
    double dx = 1.0 / 45.0;
    double dy = 0.0;
    for (int deg = 0; deg<90; deg++)
    {
        int d1 = LineHighTest(cx, cy, offsetx, offsety);
        int d2 = LineHighTest(cx, cy, 0 - offsetx, 0 - offsety);
        int d = d1;
        if (d1 > d2) d = d2;
        if (d0 > d) d0 = d;
        offsetx = offsetx + dx;
        offsety = offsety + dy;
    }
    offsetx = 1;
    offsety = 1;
    dy = 0 - dx;
    dx = 0;
    for (int deg = 0; deg<90; deg++)
    {
        int d1 = LineHighTest(cx, cy, offsetx, offsety);
        int d2 = LineHighTest(cx, cy, 0 - offsetx, 0 - offsety);
        int d = d1;
        if (d1 > d2) d = d2;
        if (d0 > d) d0 = d;
        offsetx = offsetx + dx;
        offsety = offsety + dy;
    }
    return d0;
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::CalulateFieldOpennessAllCells(int  a_pref)
{
    int dline;
    int d0 = 0;
    int minX = m_elems[a_pref]->GetMinX();
    int minY = m_elems[a_pref]->GetMinY();
    int maxX = m_elems[a_pref]->GetMaxX();
    int maxY = m_elems[a_pref]->GetMaxY();
    for (int ax = minX; ax <= maxX; ax+=10)
    {
        for (int ay = minY; ay <= maxY; ay+=10)
        {
            dline = m_maxextent; // this is the width of the landscape and will always be at least as big as the biggest return value possible
            // Get a possible point for this field
            int cx = ax;
            int cy = ay;
            if (m_land->Get(ax, ay) == a_pref)
            {
                double offsetx = -1;
                double offsety = -1;
                double dx = 1.0 / 45.0;
                double dy = 0.0;
                for (int deg = 0; deg<90; deg++)
                {
                    int d1 = LineHighTest(cx, cy, offsetx, offsety);
                    int d2 = LineHighTest(cx, cy, 0 - offsetx, 0 - offsety);
                    int d = d1;
                    if (d1 > d2) d = d2;
                    if (dline > d) dline = d; // get the minimum
                    offsetx = offsetx + dx;
                    offsety = offsety + dy;
                }
                offsetx = 1;
                offsety = 1;
                dy = 0 - dx;
                dx = 0;
                for (int deg = 0; deg<90; deg++)
                {
                    int d1 = LineHighTest(cx, cy, offsetx, offsety);
                    int d2 = LineHighTest(cx, cy, 0 - offsetx, 0 - offsety);
                    int d = d1;
                    if (d1 > d2) d = d2;
                    if (dline > d) dline = d;
                    offsetx = offsetx + dx;
                    offsety = offsety + dy;
                }
                if (dline > d0) d0 = dline; // Get the maximum. Here we might also want to do something like create statistics from the range of dline
            }
        }
    }
    return d0;
}
//-----------------------------------------------------------------------------------------------------------------------------

inline int Landscape::LineHighTest(int a_cx, int a_cy, double a_offsetx, double a_offsety)
{
    int d1=-1;
    int counter = 1;
    bool found = false;
    while (!found)
    {
        int x = (int) (a_cx + a_offsetx * counter);
        int y = (int) (a_cy + a_offsety * counter);
        if (x<1 || x >= (m_width-2) || y<1 || y >= (m_height-2)) return counter;
        TTypesOfLandscapeElement tole = this->SupplyElementType(x,y);
        if ((tole == tole_LargeRoad) || (tole == tole_SmallRoad) || (tole == tole_HedgeBank)) return counter;
        if (SupplyLEHigh(x,y))
        {
            x = (int) (a_cx + a_offsetx * (counter+1));
            y = (int) (a_cy + a_offsety * (counter+1));
            if (SupplyLEHigh(x,y)) found = true;
            d1=counter;
        }
        counter++;
    }
    return d1;
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyRodenticide(int a_x, int a_y) {
    if (cfg_rodenticide_enable.value())
    {
        double pp;
        pp = m_RodenticideManager->GetRodenticide(a_x, a_y);
        return pp;
    }
    return 0;
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyPesticide(int a_x, int a_y) {
    double pp;
#ifdef __DETAILED_PESTICIDE_FATE
    // This should be avoided if possible because it is inefficient
    pp = g_pest->SupplyPesticideS(a_x, a_y);
    pp += g_pest->SupplyPesticideP(a_x, a_y);
#else
    pp = g_pest->SupplyPesticide(a_x, a_y);
#endif
    return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::SupplyOverspray(int a_x, int a_y)
{
    return m_elems[m_land->Get(a_x, a_y)]->GetSprayedToday();
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyPesticideP(int a_x, int a_y) {
    double pp;
#ifdef __DETAILED_PESTICIDE_FATE
    pp = g_pest->SupplyPesticideP(a_x, a_y);
#else
    pp = g_pest->SupplyPesticide(a_x, a_y);
#endif
    return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyPesticideS(int a_x, int a_y) {
    double pp;
#ifdef __DETAILED_PESTICIDE_FATE
    pp = g_pest->SupplyPesticideS(a_x, a_y);
#else
    pp = g_pest->SupplyPesticide(a_x, a_y);
#endif
    return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyPesticide(int a_polyref) {
    if (!l_pest_enable_pesticide_engine.value()) return 0.0;
    double pp;
#ifdef __DETAILED_PESTICIDE_FATE
    pp = g_pest->SupplyPesticideS(a_polyref);
    pp += g_pest->SupplyPesticideP(a_polyref);
#else
    pp = g_pest->SupplyPesticide(a_polyref);
#endif
    return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyPesticideS(int a_polyref) {
    if (!l_pest_enable_pesticide_engine.value()) return 0.0;
    double pp;
#ifdef __DETAILED_PESTICIDE_FATE
    pp = g_pest->SupplyPesticideS(a_polyref);
#else
    pp = g_pest->SupplyPesticide(a_polyref);
#endif
    return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------

double Landscape::SupplyPesticideP(int a_polyref) {
    if (!l_pest_enable_pesticide_engine.value()) return 0.0;
    double pp;
#ifdef __DETAILED_PESTICIDE_FATE
    pp = g_pest->SupplyPesticideP(a_polyref);
#else
    pp = g_pest->SupplyPesticide(a_polyref);
#endif
    return pp;
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::SubtractPondLarvalFood(double a_food,int a_polyrefindex) {
    return dynamic_cast<Pond*>(m_elems[a_polyrefindex])->SubtractLarvalFood(a_food);
}
//-----------------------------------------------------------------------------------------------------------------------------

GooseFieldList* Landscape::GetGooseFields(double a_minopenness)
{
    GooseFieldList* alist = new GooseFieldList;
    GooseFieldListItem gfli;
    for (unsigned int i = 0; i < m_elems.size(); i++)
    {
        if (m_elems[i]->GetOpenness() > a_minopenness)
        {
            for (int g = gs_Pinkfoot; g < gs_foobar; g++)
            {
                gfli.grass[g] = m_elems[i]->GetGooseGrazingForage((GooseSpecies)g);
                gfli.geesesp[g] = m_elems[i]->GetGooseSpNosToday((GooseSpecies)g);
                gfli.geesespTimed[g] = m_elems[i]->GetGooseSpNosTodayTimed((GooseSpecies)g);
            }
            gfli.grain = m_elems[i]->GetBirdSeed();
            gfli.maize = m_elems[ i ]->GetBirdMaize(); 
            gfli.openness = m_elems[ i ]->GetOpenness();
            int pref = m_elems[ i ]->GetPoly();
            gfli.polyref = pref;
            gfli.geese = m_elems[i]->GetGooseNosToday();
            gfli.geeseTimed = m_elems[i]->GetGooseNosTodayTimed();
            gfli.vegtype = m_elems[i]->GetVegType();
            gfli.vegtypechr = VegtypeToString(m_elems[i]->GetVegType());
            gfli.vegheight = m_elems[i]->GetVegHeight();
            gfli.digestability = m_elems[i]->GetDigestability();
            gfli.vegphase = m_elems[i]->GetVegPhase();
            gfli.previouscrop = VegtypeToString( m_elems[ i ]->GetPreviousCrop( m_elems[ i ]->GetRotIndex() ) );
            gfli.lastsownveg = VegtypeToString( m_elems[ i ]->GetLastSownVeg() );
            alist->push_back(gfli);
        }
    }
    return alist;
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::SupplyFarmAnimalCensus(int a_farm_ref, int a_LifeStage) 
{ 
    return m_ThePopManager->FarmAnimalCensus(a_farm_ref, a_LifeStage); 
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::SupplyFarmIntensityI( int a_polyindex ) {
    return m_elems[ a_polyindex ]->GetOwner()->GetIntensity();
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::SupplyFarmIntensity( int a_x, int a_y ) {
    return m_elems[ m_land->Get( a_x, a_y ) ]->GetOwner()->GetIntensity();
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::SupplyFarmIntensity( int a_polyref ) {
  return m_elems[ m_polymapping[ a_polyref  ]]->GetOwner()->GetIntensity();
}
//-----------------------------------------------------------------------------------------------------------------------------

APoint Landscape::SupplyCentroid( int a_polyref ) {
    return m_elems[ m_polymapping[ a_polyref ] ]->GetCentroid();
}
//-----------------------------------------------------------------------------------------------------------------------------

APoint Landscape::SupplyCentroidIndex( int a_polyrefindex ) {
    return m_elems[ a_polyrefindex ]->GetCentroid();
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::SupplyPesticideDecay() 
{
    return g_pest->GetAnythingToDecay(); 
}
//-----------------------------------------------------------------------------------------------------------------------------

polylist* Landscape::SupplyLargeOpenFieldsNearXY(int a_x, int a_y, int a_range, int a_openness)
{
    polylist* p_list = new polylist;
    unsigned sz = (unsigned) m_elems.size();
    for (unsigned i=0; i<sz; i++)
    {
        if ( m_elems[i]->GetOpenness() > a_openness)
        {
            APoint pt = m_elems[i]->GetCentroid();
            int dx, dy;
            if (a_x>pt.m_x) dx = a_x-pt.m_x; else dx = pt.m_x-a_x;
            if (a_y>pt.m_y) dy = a_y-pt.m_y; else dy = pt.m_y-a_y;
            dx++; dy++; // prevents crash with maths below.
            int dist;
            if (dx>dy) dist = dx + (dy * dy) /(2 * dx); else dist = dy + (dx * dx) /(2 * dy);
            if (dist<=a_range) p_list->push_back( m_elems[i]->GetPoly());
        }
    }
    return p_list;
}
//-----------------------------------------------------------------------------------------------------------------------------

Landscape::Landscape(const char * a_configfile, const char * a_errorfile) {

    // Set up operation flags
    bool didRenumber = false;
    bool didCalcCentroids = false;
    bool didConsolidate = false;
    bool didCalcOpenness = false;
    bool didCalcOther = false;
    m_NeedCentroidCalculation = false;
    m_NeedOpennessCalculation = false;
    m_DoMissingPolygonsManipulations = false;
    // Set up globals
    g_landscape_p = this;
    for (int i = 1; i <= 2000; i++) {
        g_SpeedyDivides[i] = 1 / double(i);
    }
    int x_add[8] = { 1, 1, 0, -1, -1, -1, 0, 1 }; // W,SW,S,SE,E,NE,N,NW
    int y_add[8] = { 0, -1, -1, -1, 0, 1, 1, 1 };
    for (int i = 0; i < 8; i++) {
        m_x_add[i] = x_add[i];
        m_y_add[i] = y_add[i];
    }
    //Pinkfeet
    m_GooseIntakeRateVSVegetationHeight_PF = new Polynomial2CurveClass(cfg_P1A.value(), cfg_P1B.value(), cfg_P1C.value(), cfg_P1D.value(), cfg_P1E.value(), cfg_P1F.value(), cfg_P1G.value(), cfg_P1H.value());
    //Barnacle
    m_GooseIntakeRateVSVegetationHeight_BG = new Polynomial2CurveClass(cfg_B6A.value(), cfg_B6B.value(), cfg_B6C.value(), cfg_B6D.value(), cfg_B6E.value(), cfg_B6F.value(), cfg_B6G.value(), cfg_B6H.value());
    //Greylag
    m_GooseIntakeRateVSVegetationHeight_GL = new Polynomial2CurveClass(cfg_G6A.value(), cfg_G6B.value(), cfg_G6C.value(), cfg_G6D.value(), cfg_G6E.value(), cfg_G6F.value(), cfg_G6G.value(), cfg_G6H.value());
    if (cfg_WriteCurve.value()) {
        m_GooseIntakeRateVSVegetationHeight_GL->WriteDataFile( 10 );
        m_GooseIntakeRateVSVegetationHeight_BG->WriteDataFile( 10 );
        m_GooseIntakeRateVSVegetationHeight_PF->WriteDataFile( 10 );
    }
    
    sprintf(m_versioninfo, "%d.%d.%d :: %s", version_major, version_minor, version_revision, version_date);

    if (l_map_print_version_info.value()) {
        printf("This program uses the Landscape simulator V%s\n", m_versioninfo);
    }
    //g_cfg->DumpAllSymbolsAndExit( "allsymbols.cfg" );

    // Must come first. Used by the configurator below.
    g_msg = new MapErrorMsg(a_errorfile);
    g_msg->SetWarnLevel(WARN_ALL);

    //g_cfg->DumpAllSymbolsAndExit( "allsymbols.cfg" );

    // Configurator instantiation is automatic.
    g_cfg->ReadSymbols(a_configfile);
    //g_cfg->DumpAllSymbolsAndExit( "allsymbols.cfg" );

    // For testing.
    //g_cfg->DumpPublicSymbols( "publicsymbols.cfg", CFG_PUBLIC );
    //exit(1);

    g_date = new Calendar;
    g_weather = new Weather(l_map_weather_file.value());
    g_letype = new LE_TypeClass;
    g_crops = new CropData(l_map_cropcurves_file.value());

    m_LargestPolyNumUsed = -1;

    // Reset treatment counters.
    for (int i = 0; i < last_treatment; i++) {
        m_treatment_counts[i] = 0;
    }

    m_FarmManager = new FarmManager();


    ReadPolys2(l_map_poly_file.value());

    //this carries out the optimisation
    if (cfg_OptimisingFarms.value()){
        m_FarmManager->InitFarms();
        m_FarmManager->Save_diff_farm_types_areas();
    }


    m_land = new RasterMap(l_map_map_file.value(), this);
    m_width = m_land->MapWidth();
    m_height = m_land->MapHeight();
    m_width10 = 10 * m_width;
    m_height10 = 10 * m_height;
    if (m_width > m_height) m_maxextent = m_width; else m_maxextent = m_height;

    PolysValidate(false);

    g_pest = new Pesticide(m_land, this);
    m_toxShouldSpray = false;

    // Validate polygons, ie. ensure those reference in the
    // polygon file also shows up in the map.
    PolysValidate(false);
    PolysRemoveInvalid();

    if (cfg_MaxPondSize.value() > 0) {
        // This takes any small freshwater and converts it to a pond.
        for (unsigned int i = 0; i < m_elems.size(); i++) {
            if (m_elems[i]->GetElementType() == tole_Freshwater) {
                if (m_elems[i]->GetArea() <= cfg_MaxPondSize.value()) {
                    Pond* pond = new Pond;
                    pond->DoCopy(m_elems[i]);
                    pond->SetALMaSSEleType(g_letype->BackTranslateEleTypes(tole_Pond));
                    pond->SetElementType(tole_Pond);
                    m_elems[i] = dynamic_cast<LE*>(pond);
                }
            }
        }
    }

    PolysValidate(true);

    if (m_DoMissingPolygonsManipulations)
    {
        // Find big continuous polygons
        for (unsigned int i = 0; i < m_elems.size(); i++)
        {
            if (m_elems[i]->GetElementType() == tole_Missing)
            {
                double area = m_elems[i]->GetArea();
                int areaMinrect = (m_elems[i]->GetMaxX() - m_elems[i]->GetMinX()) * (m_elems[i]->GetMaxY() - m_elems[i]->GetMinY());
                if ((areaMinrect / area > 4) || (area < 1000))
                {
                    // Unlikely to be a field, or if so a very narrow odd one. We will assume this is a missing data issue.
                }
                else
                {
                    // Big poly with more than 25% of the included rectangle covered, must be a field of some sort.
                    // create a new wasteland and swap this in to the m_elems, then delete the old missing polygon
                    LE * wl = NewElement(tole_Wasteland);
                    wl->SetPoly(m_elems[i]->GetPoly());
                    wl->SetArea(floor(0.5 + area));
                    wl->SetSoilType(m_elems[i]->GetSoilType());
                    wl->SetUnsprayedMarginPolyRef(-1);
                    wl->SetCentroid(-1, -1);
                    wl->SetOpenness(0);
                    delete m_elems[i];
                    m_elems[i] = wl;
                }
            }
        }
        // By here all the big ones should be safely tidied away to wasteland and now we need to deal with the raster map.
        RemoveMissingValues();
        for (unsigned int i = 0; i < m_elems.size(); i++)
        {
            // Now we deal with all the little ones that were not by fields
            if (m_elems[i]->GetElementType() == tole_Missing)
            {
                LE * wl = NewElement(tole_Wasteland);
                wl->SetPoly(m_elems[i]->GetPoly());
                wl->SetArea(m_elems[i]->GetArea());
                wl->SetSoilType(m_elems[i]->GetSoilType());
                wl->SetUnsprayedMarginPolyRef(-1);
                wl->SetCentroid(-1, -1);
                wl->SetOpenness(0);
                delete m_elems[i];
                m_elems[i] = wl;
            }
        }
        PolysValidate(false);
        PolysRemoveInvalid();
        PolysValidate(true);
        g_msg->Warn("Landscape::Landscape(): Dump and normal exit to follow after resolving missing polygons. ", "");
        didCalcOther = true;
    }

    // ChangeMapMapping() also enters a valid starting
    // coordinate for the border generating farm method below.
    ChangeMapMapping();

    if ((m_LargestPolyNumUsed != ((int)m_elems.size() - 1)))
    {
        PolysRenumber();
        didRenumber = true;
    }
    // do we want to remove small polygons?
    if (l_map_removesmallpolygons.value())
    {
        cout << "In Landscape::Landscape() Small polygon removal" << endl;
        int removed = RemoveSmallPolygons();
        g_msg->Warn("Landscape::Landscape(): Dump and normal exit to follow after removing small polygons and map dump. Polygons removed:", removed);
        didCalcOther = true;
    }
    // Do we want to re-write the current files and consolidate polys?
    else if (l_map_consolidatepolys.value())
    {
        didConsolidate = true;
        ConsolidatePolys();
    }
    else if (g_map_le_borderremoval.value())
    {
        // Does not use centroids so is safe to use here
        BorderRemoval();
        CountMapSquares();
        ForceArea();
        g_msg->Warn(WARN_FILE, "Landscape::Landscape() - BorderRemoval "" map dump to follow.", "");
        didCalcOther = true;
    }


    // By the time we reach this point we need to have completed all major polygon removal and are ready to calculate centroids if they are needed. There
    // are two reasons for this - 1) that we did not have them in the original polyref file, 2) that we changed the map
    if (didConsolidate || didCalcOther || m_NeedCentroidCalculation)
    {
        PolysValidate(false);
        PolysRemoveInvalid();
        PolysValidate(true);
        ChangeMapMapping();
        PolysRenumber();
        CalculateCentroids();
        didCalcCentroids = true;
    }
    if (didConsolidate || didCalcOther || m_NeedCentroidCalculation || didCalcCentroids || l_map_calc_openness.value())
    {
        if (l_map_calc_openness.value()) CalculateOpenness(true);
        else CalculateOpenness(false);
        didCalcOpenness = true;
    }
    if (didCalcCentroids || didConsolidate || didCalcOpenness || didCalcOther || m_NeedCentroidCalculation || didRenumber || !m_FarmManager->GetIsRenumbered())
    {
        // We need to dump the map and polyrefs
        m_FarmManager->DumpFarmrefs();
        DumpMap(l_map_dump_map_file.value());
        PolysDump(l_map_dump_poly_file.value());
        g_msg->Warn(WARN_FILE, "Landscape::Landscape() ""Normal exit after dump.", "Remember to rename the new map and polyref file.");
        exit(0);
    }

    // Below here we have the more complicated map manipulations. These will need recalculation of centroids and openness after they are run.
    // However, we really do not want to get here with invalid centroids, hence forced dump and exit for manipulations up to this point.

    didCalcOther = false;
    // Add artificial hedgebanks to the hedges in the landscape,
    // if requested.
    if (l_map_art_hedgebanks.value()) {
        hb_Add();
        didCalcOther = true;
    }
    else if (g_map_le_borders.value())
    {
        cout << "Generating LE Borders around fields" << endl;
        cout << "Border chance = " << g_map_le_border_chance.value() << endl;
        cout << "Border width = " << g_map_le_borderwidth.value() << endl;
        // Generate border around each *farm* landscape element.
        cout << "Setting MaxMin Extents" << endl;
        SetPolyMaxMinExtents();
        cout << "Adding Borders" << endl;
        unsigned sz = (unsigned)m_elems.size();
        for (unsigned i = 0; i < sz; i++)
        {
            if (m_elems[i]->GetBorder() != NULL)
            {
                // Border around this element, so must be a farm field.
                // If the field is too small then ignore it
                if (m_elems[i]->GetArea() > g_map_le_borders_min_field_size.value())
                {
                    TTypesOfLandscapeElement  t = g_letype->TranslateEleTypes(g_map_le_borderstype.value());
                    BorderAdd(m_elems[i], t);
                }
            }
        }
        didCalcOther = true;
    }
    else   // Some special code to 'soften' the edges of orchards
        if (g_map_orchards_borders.value())
        {
            // Generate border around each *farm* landscape element.
            for (unsigned int i = 0; i < m_elems.size(); i++) {
                if (m_elems[i]->GetElementType() == tole_Orchard)
                {
                    TTypesOfLandscapeElement  t = g_letype->TranslateEleTypes(g_map_le_borderstype.value());
                    BorderAdd(m_elems[i], t);
                }
            }
            didCalcOther = true;
        }
        else
            // Unsprayed Margin Code....
            if (g_map_le_unsprayedmargins.value())
            {
                CountMapSquares();
                ForceArea();

                // Generate border around each *farm* landscape element.
                for (unsigned int i = 0; i < m_elems.size(); i++)
                {
                    if (m_elems[i]->GetUnsprayedMarginPolyRef() != -1)
                    {
                        // But not if the field is too small to have them (<1Ha)
                        if (m_elems[i]->GetArea() > 10000)
                        {
                            // Border around this element, so must be a farm field.
                            UnsprayedMarginAdd(m_elems[i]);
                        }
                        else   m_elems[i]->SetUnsprayedMarginPolyRef(-1);
                    }
                }
                didCalcOther = true;
            }
            else if (cfg_AddBeetleBanks.value())
            {
                cout << "Adding beetle banks now" << endl;
                AddBeetleBanks((TTypesOfLandscapeElement)cfg_BeetleBankType.value());
                didCalcOther = true;
            }

    if (l_map_dump_gfx_enable.value())
    {
        DumpMapGraphics(l_map_dump_gfx_file.value());
    }

    if (l_map_dump_enable.value() || didCalcOther)
    {
        CountMapSquares();
        ForceArea();
        PolysValidate(false);
        PolysRemoveInvalid();
        PolysValidate(true);
        ChangeMapMapping();
        PolysRenumber();
        CalculateCentroids();
        CalculateOpenness(l_map_calc_openness.value());
        m_FarmManager->DumpFarmrefs();
        cout << "Dumping map" << endl;
        DumpMap(l_map_dump_map_file.value());
        cout << "Dumping polygon refs file" << endl;
        PolysDump(l_map_dump_poly_file.value());
        g_msg->Warn(WARN_FILE, "Landscape::Landscape() ""Normal exit after dump.", "Remember to rename the new map and polyref file.");
        exit(0);
    }

    /*'''''''''''''''''' CIPE LANDSCAPE MAKER CODE HERE  //'''''''''''''''''''''
      if ( l_map_CIPEmaker_enable.value() ) {
      CIPELandscapeMaker();
      if ( l_map_dump_exit.value() ) {
      g_msg->Warn( WARN_FILE, "Landscape::Landscape(): ""Normal exit after map dump.", "" );
      exit( 0 );
      }
      }
      //'''''''''''''''''' CIPE LANDSCAPE MAKER CODE ABOVE //'''''''''''''''''''''
      */
    // Set the type of hedgebanks.
    int l_subtype = cfg_HedgeSubtypeMinimum.value();
    for (unsigned int i = 0; i < m_elems.size(); i++) {
        if (m_elems[i]->GetElementType() == tole_HedgeBank) {
            m_elems[i]->SetSubType(l_subtype);
            if (++l_subtype >= cfg_HedgeSubtypeMaximum.value())
                l_subtype = cfg_HedgeSubtypeMinimum.value();
        }
    }

    // And another to set the type of hedges
    // ***CJT*** 2003-12-02
    l_subtype = 0;
    for (unsigned int i = 0; i < m_elems.size(); i++) {
        if (m_elems[i]->GetElementType() == tole_Hedges) {
            m_elems[i]->SetSubType(l_subtype);
            if (++l_subtype >= 3)
                l_subtype = 0;
        }
    }
    // Count up the ponds and store them now we are finished with polygon handling.
    CreatePondList();

    cout << "Initiating farm management" << endl;
    m_FarmManager->InitiateManagement();
    g_date->Reset();

    /*
    if ( g_farm_test_crop.value() ) {
    TestCropManagement();
    exit( 0 );
    }
    */

    // Set up treatment flags
    // Reset internal state for the LE loop generator.
    // Compulsory!
    SupplyLEReset();
    // Get number of *all* landscape elements.
    int l_count = SupplyLECount();

    // Now loop through then.
    for (int i = 0; i < l_count; i++) {
        // Fetch next LE by its polygon reference number. Alternative
        // loop mechanism: This will return -1 at end-of-loop.
        int a_poly = SupplyLENext();

        // Skip uninteresting polygons by type, ownership,
        // phase of the moon, whatever.
        //    if ( these_are_not_the_droids_we_are_looking_for( a_poly )) {
        if (SupplyElementType(a_poly) != tole_Field)
            continue;

        // Example: Set x% of them to ignore insecticide of all types.
        if (random(100) < l_map_no_pesticide_fields.value()) {
            // Get current signal mask for polygon.
            LE_Signal l_signal = SupplyLESignal(a_poly);
            // Logical OR in/AND out the signals you are interested in.
            // The current signals are at the top of elements.h
            //l_signal |= LE_SIG_NO_INSECTICIDE | LE_SIG_NO_SYNG_INSECT | LE_SIG_NO_HERBICIDE | LE_SIG_NO_FUNGICIDE | LE_SIG_NO_GROWTH_REG;
            //l_signal |= LE_SIG_NO_INSECTICIDE | LE_SIG_NO_SYNG_INSECT | LE_SIG_NO_HERBICIDE;
            // Write the mask back out to the polygon.
            SetLESignal(a_poly, l_signal);
        }
    }

    l_vegtype_areas = (double *)malloc(sizeof(double) * (tov_Undefined + 1));

    if (l_vegtype_areas == NULL) {
        g_msg->Warn(WARN_BUG, "Landscape::Landscape(): Out of memory!", "");
        exit(1);
    }
    FILE * outf;
    if (cfg_dumpvegjan.value()) {
        outf = fopen(cfg_dumpvegjanfile.value(), "w");
        if (!outf) {
            g_msg->Warn(WARN_FILE, "Landscape::DumpMapInfoByArea(): ""Unable to create file", cfg_dumpvegjanfile.value());
            exit(1);
        }
        else
            fclose(outf);
    }

    if (cfg_dumpvegjune.value()) {
        outf = fopen(cfg_dumpvegjunefile.value(), "w");
        if (!outf) {
            g_msg->Warn(WARN_FILE, "Landscape::DumpMapInfoByArea(): ""Unable to create file", cfg_dumpvegjunefile.value());
            exit(1);
        }
        else
            fclose(outf);
    }

    // Dump veg information if necessary
    if (l_map_dump_veg_enable.value()) {
        FILE * f;
        f = fopen("VegDump.txt", "w");
        if (!f) {
            g_msg->Warn(WARN_BUG, "Landscape::Landscape(): VegDump.txt could not be created", "");
            exit(1);
        }
        fprintf(f, "Year\tDay\tHeight\tBiomass\tGrazed\tDensity\tCover\tWeedBiomass\ttovNum\tInsectBiomass\tLATotal\tLAGreen\tDigestability\tGreenBiomass\tDeadBiomass\tGooseGrazing\tSpilledGrain\nn");
        fclose(f);
    }
    if (l_map_dump_event_enable.value()) {
        FILE * f;
        f = fopen("EventDump.txt", "w");
        if (!f) {
            g_msg->Warn(WARN_BUG, "Landscape::Landscape(): EventDump.txt could not be created", "");
            exit(1);
        }
        fclose(f);
    }

    if (!cfg_OptimiseBedriftsmodelCrops.value()){
        m_FarmManager->FindNeighbours();
    }

    if (cfg_DumpFarmAreas.value()){
        m_FarmManager->DumpFarmAreas();
    }

    // If we are testing a pesticide then set the enum attribute
    m_PesticideType = (TTypesOfPesticide)cfg_pesticidetesttype.value();

    if (cfg_rodenticide_enable.value())
    {
        m_RodenticideManager = new RodenticideManager("BaitLocations_input.txt", this);
        m_RodenticidePreds = new RodenticidePredators_Population_Manager(this);
    }

    // Run a year to remove any start up effects
    cout << "Running intial start-up year" << endl;
    for (unsigned int i = 0; i < 365; i++)  Tick();
    //switch the rotation after running the hidden year (only for optimising farms), AM, 030713
    if (cfg_OptimisingFarms.value()) { m_FarmManager->Switch_rotation(); }
    // Write ASCII file:
    if (l_map_write_ascii.value()) {
        int x = l_map_ascii_utm_x.value();
        int y = l_map_ascii_utm_y.value();
        GISASCII_Output("AsciiLandscape.txt", x, y);
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::ConsolidatePolys()
{
        const int TypesToReplace = 18;
        TTypesOfLandscapeElement replaceList[TypesToReplace] = { tole_River, tole_RiversidePlants, tole_RiversideTrees, tole_StoneWall, tole_BareRock, tole_BuiltUpWithParkland,
            tole_Carpark, tole_Churchyard, tole_Coast, tole_Garden, tole_HeritageSite,
            tole_IndividualTree, tole_PlantNursery, tole_Saltwater, tole_SandDune, tole_UrbanNoVeg };
        int foundList[TypesToReplace];
        cout << "Consolidating polygons with no special behaviour" << endl;
        for (int i = 0; i < TypesToReplace; i++) foundList[i] = -1;
        int mapwidth = m_land->MapWidth();
        int mapheight = m_land->MapHeight();
        for (int x = 0; x < mapwidth; x++) 
        {
            for (int y = 0; y < mapheight; y++) 
            {
                int ele = m_land->Get(x, y);
                TTypesOfLandscapeElement tole = m_elems[m_polymapping[ele]]->GetElementType();
                for (int t = 0; t < TypesToReplace; t++) 
                {
                    if (tole == replaceList[t])
                    {
                        // Must do something with this cell
                        if (foundList[t] == -1) foundList[t] = ele;
                        else
                        {
                            // Need to replace this cell
                            m_land->Put(x, y, foundList[t]);
                        }
                    }
                }
            }
        }
    // At this point there should be many polygons that are not in the map. So we need to run the valid test.
    g_msg->Warn(WARN_FILE, "Landscape::ConsolidatePolys() - ""Consolidate map dump.", "");
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::DumpMap( const char * a_filename ) 
{
    int * l_map = m_land->GetMagicP(0, 0); // Hmmm - this is a nasty way round the class data protection. Gets a pointer direct to m_map in rastermap.
    /*  FILE * l_file;
  l_file = fopen(a_filename, "wb" );
  if ( !l_file ) {
    g_msg->Warn( WARN_FILE, "Landscape::DumpMap(): Unable to open file", a_filename );
    exit( 0 );
  }

  char * l_id = m_land->GetID();


      fwrite( l_id, 1, 12, l_file );
  fwrite( & m_width, 1, sizeof( int ), l_file );
  if (cfg_rectangularmaps_on.value() )
  {
      fwrite( & m_height, 1, sizeof( int ), l_file );
  }
  for ( int i = 0; i < m_width * m_height; i++ ) 
  {
        LE* le = m_elems[m_polymapping[l_map[i]]];
        int l_poly = le->GetPoly();
        fwrite( & l_poly, 1, sizeof( int ), l_file );
  }
  fclose( l_file );
  */
    ofstream OFile( a_filename, ios::binary);
  char id[12] = { "LSB2_Format" };
  OFile.write(id, 12);
  OFile.write((char*)&m_width, sizeof (int));
  OFile.write((char*)&m_height, sizeof (int));
  OFile.write((char*)l_map, m_width*m_height*sizeof (int));
  OFile.close();
}
//-----------------------------------------------------------------------------------------------------------------------------


Landscape::~Landscape( void ) {
  if ( l_map_dump_treatcounts_enable.value() ) {
    DumpTreatCounters( l_map_dump_treatcounts_file.value() );
  }

  for ( unsigned int i = 0; i < m_elems.size(); i++ )
    delete m_elems[ i ];

  free( m_polymapping );
  free( l_vegtype_areas );
  delete m_land;
  //delete g_rotation;
  delete g_crops;
  delete g_letype;
  delete g_weather;
  delete g_date;
  delete g_pest;
  if (cfg_rodenticide_enable.value())
  {
      delete m_RodenticideManager;
      delete m_RodenticidePreds;
  }
  delete m_FarmManager;
  delete g_msg; // Must be last.
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::SimulationClosingActions()
{
    if(cfg_OptimisingFarms.value()){
        this->m_FarmManager->PrintFinalResults();
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::Tick( void ) {
    g_date->Tick();
    g_weather->Tick();

    // Remember todays LAItotal for veg elements
    for (unsigned int i = 0; i < m_elems.size(); i++) {
        m_elems[ i ]->StoreLAItotal();
    }

    // Update the growth curve phases if needed.
    if (g_date->JanFirst()) {
        if (g_rand_uni() > 0.5) {
            m_FarmManager->SetSpilledGrain( true );
        }
        else m_FarmManager->SetSpilledGrain( false );
        for (unsigned int i = 0; i < m_elems.size(); i++) {
            m_elems[ i ]->SetGrowthPhase( janfirst );
        }
    }
    else if (g_date->MarchFirst()) {
        for (unsigned int i = 0; i < m_elems.size(); i++) {
            m_elems[ i ]->SetGrowthPhase( marchfirst );
        }
        // Check and see if the pesticide engine flag should be set
        if ((SupplyYearNumber() >= cfg_productapplicstartyear.value()) && (SupplyYearNumber() <= cfg_productapplicendyear.value()))
            m_toxShouldSpray = true;
        else m_toxShouldSpray = false;
    }


    // Grow the green stuff and let the bugs have some too.
    for (unsigned int i = 0; i < m_elems.size(); i++) {
        m_elems[ i ]->Tick();
        m_elems[ i ]->DoDevelopment();
    }

    //call ActualProfit function, then save this year's crop prices - for the next year's expected results calculation, then change price
    if (cfg_OptimisingFarms.value()) {
        if (!cfg_OptimiseBedriftsmodelCrops.value()) { //do this just in almass crop mode

            if (g_date->DayInYear() == g_date->DayInYear( 1, 03 ) && g_date->GetYearNumber()>1) { //March 1st - do the accounting; 1 corresponds to 0 in the visual version year counter.

                m_FarmManager->ActualProfit();

                //print the m_cropTotals
                ofstream ofile2( "CropDistribution.txt", ios::app );
                // the distibution of crops from last year! 
                int last_year = g_date->GetYearNumber() - 1;
                ofile2 << last_year << '\t';
                for (int c = 0; c < m_FarmManager->pm_data->Get_noCrops(); c++) ofile2 << m_FarmManager->Get_cropTotals( c ) << '\t';
                ofile2 << endl;
                //restart the m_cropTotals:
                for (int i = 0; i < (int)m_FarmManager->Get_cropTotals_size(); i++) {
                    m_FarmManager->Set_cropTotals( i, 0 );
                }

                if (g_date->GetYearNumber() > 7) { //it's been 6 years since the first accounting was carried out (5 would be enough, but exclude one year to avoid start up effects)
                    m_FarmManager->Save_last_years_crops_for_farms(); //saves m_rotational_crops in m_rotational_crops_visible
                    m_FarmManager->ChooseDecisionMode_for_farms();
                }

            }

            //save the last year crop price; 
            if (g_date->GetYearNumber() > 1 && g_date->DayInYear() == g_date->DayInYear( 1, 03 )) {
                for (int f = 0; f < (toof_Foobar * tos_Foobar); f++) { //need to save the price for all farm types and soil types, so just run through the whole vector;
                    int no_crops = m_FarmManager->pm_data->Get_cropTypes_almass_size();
                    for (int c = 0; c < no_crops; c++) {
                        double old_price = m_FarmManager->pm_data->Get_sellingPrice( tov_Undefined*f + m_FarmManager->pm_data->Get_cropTypes_almass( c ) ); //get the current price
                        m_FarmManager->pm_data->Set_sellingPrice_lastyr( old_price, tov_Undefined*f + m_FarmManager->pm_data->Get_cropTypes_almass( c ) ); //save it
                    }
                }
            }
            //-----


            if (cfg_MaizeEnergy.value()) { //modify energy maize price

                if (g_date->DayInYear() == g_date->DayInYear( 1, 03 ) && g_date->GetYearNumber() != 0) { //print the price once a year
                    ofstream ofileEM( "Maize_energy_price_yearly.txt", ios::app );
                    int year = g_date->GetYearNumber();
                    ofileEM << year << '\t';
                    ofileEM << m_FarmManager->pm_data->Get_sellingPrice( tov_Maize );
                    ofileEM << endl;
                    ofileEM.close();
                }

                //changing price according to the input file EM_price.txt
                if (g_date->DayInYear() == g_date->DayInYear( 1, 03 ) && g_date->GetYearNumber() > 0) {//exclude the hidden year
                    for (int f = 0; f < (toof_Foobar * tos_Foobar); f++) { //need to change the price for all farm types and soil types
                        double new_price = m_FarmManager->pm_data->Get_emaize_price( g_date->GetYearNumber() );
                        m_FarmManager->pm_data->Set_sellingPrice( new_price, tov_Undefined*f + tov_Maize );
                    }
                }
            }// if energy maize

        }//almass crop mode
    }//optimising farms

    /*  Testing removal - may cause issues with increasing or decreasing permanent vegetation
    if ( g_date->DayInYear() == g_date->DayInYear( 1, 11 ) ) {
    // Set all elements to smooth curve transition mode at November 1st.
    for ( unsigned int i = 0; i < m_elems.size(); i++ ) {
    m_elems[ i ]->ForceGrowthInitialize();
    }
    //save the day degrees - used for crops that are not harvested instead of biomass to determine yield; done on Nov 1st
    //m_FarmManager->SetDD (SupplyTempPeriod( g_date->Date(), 245)); // from Nov 1 - 245 days back till March 1;
    //DegreesDump();
    }
    */

    // Put the farmers to work.
    m_FarmManager->FarmManagement();


    // Update pesticide information.
    g_pest->Tick();

    // Update rodenticide information if we are using this
    if (cfg_rodenticide_enable.value()) {
        m_RodenticideManager->Tick();
        m_RodenticidePreds->Tick();
    }

    // Dump event information if necessary
    if (l_map_dump_veg_enable.value()) VegDump( l_map_dump_veg_x.value(), l_map_dump_veg_y.value() );
    if (l_map_dump_event_enable.value())
        //EventDump( l_map_dump_event_x1.value(), l_map_dump_event_y1.value(), l_map_dump_event_x2.value(), l_map_dump_event_y2.value() );
        EventDumpPesticides( l_map_dump_event_x1.value(), l_map_dump_event_y1.value() );
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::SetPolyMaxMinExtents()
{
  // All polygon manipulation is settled now we need to give the polygons some information about themselves
  // This takes a little time but save time later one
    cout << "Setting max min polygon extents" << endl;
    int mwidth = m_land->MapWidth();
    int mheight = m_land->MapHeight();
    for ( int x = 0; x < mwidth; x++ ) {
        for ( int y = 0; y < mheight; y++ ) {
            int polyindex = m_land->Get( x, y );
            // Mark that we have seen this polygon.
            unsigned int ele_ref= polyindex;
            if (m_elems[m_polymapping[ele_ref]]->GetMaxX() < x) m_elems[m_polymapping[ele_ref]]->SetMaxX(x);
            if (m_elems[m_polymapping[ele_ref]]->GetMaxY() < y) m_elems[m_polymapping[ele_ref]]->SetMaxY(y);
            if (m_elems[m_polymapping[ele_ref]]->GetMinX() > x) m_elems[m_polymapping[ele_ref]]->SetMinX(x);
            if (m_elems[m_polymapping[ele_ref]]->GetMinY() > y) m_elems[m_polymapping[ele_ref]]->SetMinY(y);
            m_elems[m_polymapping[ele_ref]]->SetMapValid(true);
        }
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::VegDump( int x, int y ) {
  FILE * vfile=fopen("VegDump.txt", "a" );
  if (!vfile) {
    g_msg->Warn( WARN_FILE, "Landscape::VegDump(): Unable to open file", "VegDump.txt" );
    exit( 1 );
  }
  int year = SupplyYearNumber();
  int day = SupplyDayInYear();
  double hei = SupplyVegHeight( x, y );
  double bio = SupplyVegBiomass( x, y );
  double cover = SupplyVegCover( x, y );
  double density = bio / ( hei + 1 );
  double weeds = SupplyWeedBiomass( x, y );
  double insects = SupplyInsects( x, y );
  double LATotal = SupplyLATotal(x, y);
  double LAGreen = SupplyLAGreen(x, y);
  double digest = SupplyVegDigestability(x, y);
  double GreenBiomass = SupplyGreenBiomass(x,y);
  double DeadBiomass = SupplyDeadBiomass(x,y);
  int grazed = SupplyGrazingPressure(x, y);
  double ggraze = GetActualGooseGrazingForage(m_land->Get(x, y), gs_Pinkfoot);
  int VegType = BackTranslateVegTypes(SupplyVegType(x, y));
  double grain = SupplyBirdSeedForage(x, y);

  fprintf( vfile, "%d\t%d\t%g\t%g\t%d\t%g\t%g\t%g\t%i\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", year, day, hei, bio, grazed, density, cover, weeds, VegType, insects, LATotal, LAGreen, digest, GreenBiomass, DeadBiomass, ggraze, grain );
  /*
  y += 100;
  x += 100;
  hei = SupplyVegHeight( x, y );
  bio = SupplyVegBiomass( x, y );
  cover = SupplyVegCover( x, y );
  density = bio / ( hei + 1 );
  weeds = SupplyWeedBiomass( x, y );
  VegType = SupplyVegType( x, y );
  insects = SupplyInsects( x, y );
  digest = SupplyVegDigestability(x,y);
  fprintf( vfile, "%g\t%g\t%g\t%g\t%g\t%i\t%g\t%g\n", hei, bio, density, cover, weeds, VegType, insects, digest );
  */
  fclose( vfile );
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::DegreesDump(){

    ofstream ofile ("Daydegrees.txt", ios::app);    
    //print degrees
    ofile << m_FarmManager->GetDD();
    ofile  << endl;
    ofile.close();
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::EventDump( int x1, int y1, int x2, int y2 ) {
  FILE * vfile=fopen("EventDump.txt", "a" );
  if (!vfile) {
    g_msg->Warn( WARN_FILE, "Landscape::EventDump(): Unable to open file", "EventDump.txt" );
    exit( 1 );
  }
  FarmToDo event;
  int i = 0;
  int day = SupplyDayInYear();
  fprintf( vfile, "%d: ", day );
  while ( ( event = ( FarmToDo )SupplyLastTreatment( x1, y1, & i ) ) != sleep_all_day ) {
    fprintf( vfile, "%d ", event );
  }
  i = 0;
  fprintf( vfile, " - " );
  while ( ( event = ( FarmToDo )SupplyLastTreatment( x2, y2, & i ) ) != sleep_all_day ) {
    fprintf( vfile, "%d ", event );
  }
  fprintf( vfile, "\n" );
  fclose( vfile );
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::EventDumpPesticides( int x1, int y1 ) {
  FILE * vfile=fopen("EventDump.txt", "a" );
  if (!vfile) {
    g_msg->Warn( WARN_FILE, "Landscape::EventDump(): Unable to open file", "EventDump.txt" );
    exit( 1 );
  }
  FarmToDo a_event;
  int i = 0;
  int day = this->SupplyGlobalDate();
  int herb = 0;
  int fung = 0;
  int ins = 0;
  while ( ( a_event = ( FarmToDo )SupplyLastTreatment( x1, y1, & i ) ) != sleep_all_day ) {
      if (a_event == herbicide_treat )
      {
              herb++;
      }
      else if (a_event == fungicide_treat ) fung++;
      else if (a_event == insecticide_treat) ins++;
  }
  if (herb+fung+ins >0 ) fprintf( vfile, "%d\t%d\t%d\t%d\n", day, herb, fung, ins );
  i = 0;
  fclose( vfile );
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::PolysValidate( bool a_exit_on_invalid ) {
   // First loop just sets the MapValid as false (and checks for a major screw-up if this elemenent does not exist even in the list
  for ( unsigned int i = 0; i < m_elems.size(); i++ ) {
    m_elems[ i ]->SetMapValid( false );
    if ( m_polymapping[ m_elems[ i ]->GetPoly() ] == -1 ) {
      char l_err[ 20 ];
      sprintf( l_err, "%d", m_elems[ i ]->GetPoly() );
      g_msg->Warn( WARN_FILE, "Landscape::PolysValidate(): Invalid polymapping ", l_err );
      exit( 1 );
    }
  }
  // Now go through the whole map and for each polygon found set MapValid as true.
  SetPolyMaxMinExtents();

  if ( a_exit_on_invalid ) {
    for ( unsigned int i = 0; i < m_elems.size(); i++ ) {
      if ( !m_elems[ i ]->GetMapValid() ) {
        char l_err[ 20 ];
        sprintf( l_err, "%d", m_elems[ i ]->GetPoly() );
        g_msg->Warn( WARN_FILE, "Landscape::PolysValidate(): Invalid polygon ", l_err );
        exit( 0 );
      }
    }
  }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::PolysRemoveInvalid( void ) 
{
  vector < LE * > l_temp;
  cout << "Tidying up the polygon map in PolysRemoveInvalid" << endl;
  unsigned int sz= (int) m_elems.size();
  for ( unsigned int i = 0; i < sz; i++ ) {
    if ( m_elems[ i ]->GetMapValid() ) {
      unsigned int j =  (int) l_temp.size();
      l_temp.resize( j + 1 );
      l_temp[ j ] = m_elems[ i ];
    } else {
        // cout << "Deleted m_elems index:" << m_polymapping[ m_elems[ i ]->GetPoly() ] << " Polynumber :" << m_elems[ i ]->GetPoly() << BackTranslateEleTypes(m_elems[ i ]->GetElementType()) << endl;
        m_polymapping[ m_elems[ i ]->GetPoly() ] = -1;
      delete m_elems[ i ];
    }
  }

  for ( unsigned int i = 0; i < l_temp.size(); i++ ) {
    m_elems[ i ] = l_temp[ i ];
  }
  m_elems.resize( l_temp.size() );
  RebuildPolyMapping();
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::PolysDump(const char * a_filename)
{
    ofstream outf(a_filename, ios::out);
    int l_num_polys = 0;

    if (!outf.is_open()) {
        g_msg->Warn(WARN_FILE, "Landscape::PolysDump(): Unable to open file", a_filename);
        exit(1);
    }

    // Count up number if active polygons in our list.
    unsigned sz = (unsigned)m_elems.size();
    for (unsigned int i = 0; i < sz; i++) {
        if (m_elems[i]->GetMapValid())
            l_num_polys++;
    }

    outf << l_num_polys << endl;
    outf << "PolyType" << '\t' << "PolyRefNum" << '\t' << "Area" << '\t' << "FarmRef" << '\t' << "UnSprayedMarginRef" << '\t' << "SoilType" << '\t' << "Openness" << '\t' << "CentroidX" << '\t' << "CentroidY" << endl;

    /*
    if ( l_map_renumberpolys.value() )
    {
    RenumberPolys( true );
    }
    */
    // Now we can output the file
    for (unsigned int i = 0; i < m_elems.size(); i++)
    {
        if (m_elems[i]->GetMapValid())
        {
            outf << m_elems[i]->GetALMaSSEleType() << '\t' << m_elems[i]->GetPoly() << '\t' << m_elems[i]->GetArea() << '\t' <<
                m_elems[i]->GetOwnerFile() << '\t' << m_elems[i]->GetUnsprayedMarginPolyRef() << '\t' << m_elems[i]->GetSoilType() << '\t' << m_elems[i]->GetOpenness()
                << '\t' << m_elems[i]->GetCentroidX() << '\t' << m_elems[i]->GetCentroidY() << endl;
        }
    }
    outf.close();
}
//-----------------------------------------------------------------------------------------------------------------------------

/*
void Landscape::RenumberPolys( bool a_checkvalid)
{
      // First need a second list of polyrefs
      vector < int > oldpolys;
      oldpolys.resize( m_elems.size() );
      if (a_checkvalid)
      {
          for ( unsigned int i = 0; i < m_elems.size(); i++ )
          {
              if ( m_elems[ i ]->GetMapValid() )
              {
                  oldpolys[ i ] = m_elems[ i ]->GetPoly();
                  // Need to reset the poly number
                  m_elems[ i ]->SetPoly( i );
              }
          }
      }
      else
      {
          for ( unsigned int i = 0; i < m_elems.size(); i++ )
          {
              oldpolys[ i ] = m_elems[ i ]->GetPoly();
              // Need to reset the poly number
              m_elems[ i ]->SetPoly( i );
          }
      }

      // Now need to go through and set the unsprayed margin poly refs
      for ( unsigned int j = 0; j < m_elems.size(); j++ )
      {
          int um = m_elems[ j ]->GetUnsprayedMarginPolyRef();
          if ( um != -1 )
          {
              unsigned ff =  (int) oldpolys.size();
              for ( unsigned int u = 0; u < ff; u++ )
              {
                  if ( oldpolys[ u ] == um )
                  {
                      m_elems[ j ]->SetUnsprayedMarginPolyRef( u );
                      break;
                  }
              }
          }
      }
}
//-----------------------------------------------------------------------------------------------------------------------------
*/

void Landscape::ReadPolys2(const char * a_polyfile)
{
    // Need to figure out if the farms are already renumbered - as the Farm manager
    bool farmsrenum = m_FarmManager->GetIsRenumbered();

    int NoPolygons;
    string rubbish = ""; //put here the names of the parameters;
    ifstream ifile(a_polyfile);
    if (!ifile.is_open()) {
        g_msg->Warn(WARN_FILE, "Landscape::ReadPolys(): Unable to open file", a_polyfile);
        std::exit(1);
    }
    char error_num[20];
    // First read the number of polygons 
    ifile >> NoPolygons;
    m_elems.resize(NoPolygons);
    for (int i = 0; i < 9; i++){ ifile >> rubbish; }
    int np = NoPolygons;
    if (NoPolygons < 10000) np = 10000; // this is just to cope with old maps with < 10000 polygons that do not follow the rules for new maps
    m_polymapping = (int *)malloc(sizeof(int)* (np * 2));
    // Test if we have the memory to do this
    if (m_polymapping == NULL) 
    {
        g_msg->Warn(WARN_FILE, "Landscape::ReadPolys(): Out of memory!", "");
        std::exit(1);
    }

    // Set all mappings to unused.
    for (int i = 0; i < np * 2; i++) {
        m_polymapping[i] = -1;
    }

    int ElemIndex = 0;

    for (int x = 0; x < NoPolygons; x++)
    {
        int PolyNum, Owner, PolyType, RealPolyType, URef, SoilType, openness, Centroid_x, Centroid_y;
        TTypesOfLandscapeElement Type;
        float Area;
        ifile >> PolyType >> PolyNum >> Area >> Owner >> URef >> SoilType >> openness >> Centroid_x >> Centroid_y;
        // Here we make some tests to check input validity
        if ((SoilType > 16) || (PolyNum<0))
        {
            std::sprintf(error_num, "%d", NoPolygons);
                g_msg->Warn(WARN_FILE, "Landscape::ReadPolys(): Polygon file empty before "
                    "reading number of specified polygons (old polygon file format?):", error_num);
                std::exit(1);
        }

        // Owner is the farm number or -1. If farms have not been renumbered then this needs mapped to the index in the list of farms.
        // Because we create this sequentially in Farm.cpp we have a constant index.
        if ((-1 != Owner) && !farmsrenum)
        {
            // Need to replace the Owner with the new renumbered ref.
            Owner = m_FarmManager->GetRenumberedFarmRef(Owner);
        }

        RealPolyType = PolyType;
        if (PolyType == 150) 
        {
            PolyType = l_map_chameleon_replace_num.value();
        }

        Type = g_letype->TranslateEleTypes(PolyType);
        if (Type == tole_Missing)
        {
            m_DoMissingPolygonsManipulations = true;
        }
        if (-1 == m_polymapping[PolyNum])
        {
            // First time we have encountered this polygon number.
            // Borders are not mapped in this list.
            m_polymapping[PolyNum] = ElemIndex;
            LE * newland = NewElement(Type);
            m_elems[ElemIndex++] = newland;
            newland->SetPoly(PolyNum);
            newland->SetArea(floor(0.5 + Area));
            newland->SetALMaSSEleType(RealPolyType);
            newland->SetSoilType(SoilType);
            newland->SetUnsprayedMarginPolyRef(URef);
            newland->SetCentroid(Centroid_x,Centroid_y);
            newland->SetOpenness(openness);
            // Just for fun, or maybe because we might need it later, remember the actual largest polynum used
            if (PolyNum>m_LargestPolyNumUsed) m_LargestPolyNumUsed = PolyNum;
            // Now set any centroid or openness recalcuation flags
            if ((Centroid_x < 0) || (Centroid_y < 0)) m_NeedCentroidCalculation= true;
            if (openness < 0) m_NeedOpennessCalculation = true;
            // Two types of possible errors: Landscape element that is a field,
            // but doesn't belong to a farm, or a farm element not of type field.
            // Check for both cases.
            if (-1 == Owner && (Type == tole_Field || Type == tole_YoungForest || Type == tole_Orchard || Type == tole_PermPastureTussocky || Type == tole_PermPasture || Type == tole_PermanentSetaside || Type == tole_PermPastureLowYield || Type == tole_WoodyEnergyCrop || Type == tole_Vildtager || Type == tole_PlantNursery)) {
                // No owner but field polygon.
                sprintf(error_num, "%d", PolyNum);
                g_msg->Warn(WARN_FILE, "Landscape::ReadPolys(): Farm polygon does not belong to a farm:", error_num);
                exit(1);
            }
            if (-1 != Owner && Type != tole_Field && Type != tole_YoungForest && Type != tole_Orchard && Type != tole_PermPastureTussocky && Type != tole_PermPasture && Type != tole_PermanentSetaside && Type != tole_PermPastureLowYield && Type != tole_WoodyEnergyCrop && Type != tole_Vildtager && Type != tole_PlantNursery) {
                // An owner but not field elements.
                sprintf(error_num, "%d", PolyNum);
                g_msg->Warn(WARN_FILE, "Landscape::ReadPolys(): Farm polygon does not have element type tole_Field:", error_num);
                exit(1);
            }

            if (-1 != Owner)
            {
                m_FarmManager->ConnectFarm(Owner);
                m_FarmManager->AddField(Owner, newland, Owner);
                if (g_map_le_borders.value())
                {
                    if (random(100) < g_map_le_border_chance.value())
                    {
                        // This is a farm element, so signal adding a border.
                        newland->SetBorder((LE *)1);
                    }
                }
                // Code to generate unsprayed margins....
                if (newland->GetElementType() == tole_Field)
                {
                    if (g_map_le_unsprayedmargins.value())
                    {
                        if (random(100) < g_map_le_unsprayedmargins_chance.value())
                        {
                            // This is a farm field, so signal adding a margin
                            newland->SetUnsprayedMarginPolyRef(1);
                        }
                    }
                }
                // ..to here
            }
        }
    else {
            sprintf(error_num, "%d", PolyNum);
            g_msg->Warn(WARN_FILE, "Landscape::ReadPolys(): Duplicate polygon in file", error_num);
            exit(1);
        }
    }
    ifile.close();
    hb_first_free_poly_num = m_elems[NoPolygons - 1]->GetPoly() + 1;
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::CountMapSquares( void ) {
    int mapwidth = m_land->MapWidth();
    int mapheight = m_land->MapHeight();
    for ( unsigned int i = 0; i < m_elems.size(); i++ ) {
        m_elems[i]->SetArea(0);
        m_elems[ i ]->m_squares_in_map=0;
    }

    for ( int x = 0; x < mapwidth; x++ ) {
        for ( int y = 0; y < mapheight; y++ ) {
            int l_ele = m_land->Get( x, y );
            m_elems[ l_ele ]->m_squares_in_map++;
        }
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::PolysRenumber(void)
{
    cout << "In Landscape::Landscape() Polygon renumber." << endl;
    for (unsigned int i = 0; i < m_elems.size(); i++)
    {
        // Need to reset the poly number
        int index = m_elems[i]->GetMapIndex(); // This is the number currently in the map matrix
        m_elems[i]->SetPoly(index); // The map index and the polygon number are now one and the same
        m_polymapping[index] = i; // The polymapping is now linked via index to the m_elems index (i)
    }
    m_LargestPolyNumUsed = (int) m_elems.size()-1;
    g_msg->Warn(WARN_FILE, "Landscape::Landscape() ""Map to be dumped due to polygon renumber", "");
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::ForceArea( void ) {
  int l_area_sum = 0;

  for ( unsigned int i = 0; i < m_elems.size(); i++ ) {
    m_elems[ i ]->SetArea( ( double )m_elems[ i ]->m_squares_in_map );
    if ( m_elems[ i ]->m_squares_in_map > 0 ) {
      m_elems[ i ]->SetMapValid( true );
      l_area_sum += m_elems[ i ]->m_squares_in_map;
    }
  }

  if ( l_area_sum != m_width * m_height ) {
    g_msg->Warn( WARN_BUG, "Landscape::ForceArea(): Polygon areas doesn't"" sum up to map area!", "" );
    exit( 1 );
  }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::RemoveMissingValues(void)
{
    bool found; // the flag for something left to work with
    int mapwidth = m_land->MapWidth();
    int mapheight = m_land->MapHeight();
    int counter = 0;
    do
    {
        found = false; counter++;
        for (int x = 1; x < mapwidth-1; x++)
        {
            for (int y = 1; y < mapheight-1; y++)
            {
                int apoly = m_land->Get(x,y);
                if (m_elems[m_polymapping[apoly]]->GetElementType() == tole_Missing)
                {
                    m_land->MissingCellReplace(x, y, true);
                }
            }
        }
        counter++;
    } while (counter<50);
    // Now we only have the edges to deal with 
    for (int x = 0; x < mapwidth; x++)
    {
        for (int y = 0; y < mapheight; y++)
        {
            int apoly = m_land->Get(x, y);
            if (m_elems[m_polymapping[apoly]]->GetElementType() == tole_Missing)
            {
                found = true;
                counter++;
                m_land->MissingCellReplaceWrap(x, y, true);
            }
        }
    }
    // Now we the ones that are not next to fields to deal with 
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::ChangeMapMapping( void ) {
    cout << "In Change Map Mapping" << endl;
    int mapwidth = m_land->MapWidth();
    int mapheight = m_land->MapHeight();
    int pest_map_width = mapwidth >> PEST_GRIDSIZE_POW2;
    if ( mapwidth & ( PEST_GRIDSIZE - 1 ) ) pest_map_width++;
    int oldindex = -1;
    for ( int x = 0; x < mapwidth; x++ ) 
    {
        for ( int y = 0; y < mapheight; y++ ) 
        {
            int polynum = m_land->Get( x, y ); // the polyref e.g. = 1, m_polymapping[ polynum ] = 0
            m_elems[ m_polymapping[ polynum ]]->SetMapIndex( m_polymapping[ polynum ] ); // Here we set index in the map to the index in elements, i.e. 0
            m_elems[ m_polymapping[ polynum ]]->SetMapValid( true );
            // Do the translation.
            m_land->Put( x, y, m_polymapping[ polynum ] ); // and now we write this to the map, i.e. 0
            // This coordinate is now valid. Throw these coordinates into
            // the associated landscape element.
            int index = m_polymapping[ SupplyPolyRef( x, y ) ];
            if ( index != oldindex ) 
            {
                m_elems[ index ]->SetValidXY( x, y );
                int l_x = x >> PEST_GRIDSIZE_POW2;
                int l_y = y >> PEST_GRIDSIZE_POW2;
                int pref = l_y * pest_map_width + l_x;        
                m_elems[ index ]->SetPesticideCell( pref );
                oldindex = index;
            }
        }
    }
    RebuildPolyMapping(); 
/*
// Check that all of the polygons are mentioned in the map.
    if ( l_map_check_polygon_xref.value() ) 
    {
        for ( unsigned int i = 0; i < m_elems.size(); i++ ) 
        {
            if ( !m_elems[ i ]->GetMapValid() ) {
                char poly[ 20 ];
                sprintf( poly, "%d", m_elems[ i ]->GetPoly() );
                g_msg->Warn( WARN_FILE, "Landscape::ChangeMapMapping(): ""Polygon number referenced but not in map file: ", poly );
                exit( 1 );
            }
        }
    }
*/
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::RemoveSmallPolygons()
{
    // To be sure we first count the map squares
    CountMapSquares();
    // First find polygons that don't match other rules and are single cells - remove these
    int removed = 0;
    for (unsigned int i = 0; i < m_elems.size(); i++)
    {
        if (m_elems[i]->GetArea() == 1)
        {
            TTypesOfLandscapeElement tole = m_elems[i]->GetElementType();
            switch (tole)
            {
                case tole_FieldBoundary:
                case tole_HedgeBank:
                case tole_Hedges:
                case tole_IndividualTree:
                case tole_MetalledPath:
                case tole_River:
                case tole_RiversidePlants:
                case tole_RiversideTrees:
                case tole_RoadsideSlope:
                case tole_RoadsideVerge:
                case tole_SmallRoad:
                case tole_StoneWall:
                case tole_Stream:
                case tole_Track:
                case tole_UnsprayedFieldMargin:
                    // These could be part of a bigger polygon, so leave them alone for now.
                    break;
                default:
                    // Get rid of this one.
                    APoint pt = m_elems[i]->GetCentroid();
                    m_land->CellReplacementNeighbour(pt.m_x, pt.m_y, m_elems[i]->GetPoly());
                    removed++;
            }
        }
    }
    return removed;
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::BorderRemoval() {
    // This does not need to be efficient, just do the job
    for (int x=1; x<(m_width-1); x++)
        for (int y=1; y<(m_height-1); y++)
        {
            TTypesOfLandscapeElement tole = SupplyElementType(x,y);
            if ((tole==tole_FieldBoundary) || (tole==tole_HedgeBank) || (tole==tole_Hedges))
            {
                if ( SupplyElementType(x-1,y-1) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x-1,y-1);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x-1,y) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x-1,y);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x-1,y+1) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x-1,y+1);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x,y-1) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x,y-1);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x,y+1) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x,y+1);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x+1,y-1) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x+1,y-1);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x+1,y) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x+1,y);
                    m_land->Put( x, y, fieldindex );

                }
                else
                if ( SupplyElementType(x+1,y+1) == tole_Field)
                {
                    // Set the x,y location to be this field
                    int fieldindex = SupplyPolyRefIndex(x+1,y+1);
                    m_land->Put( x, y, fieldindex );

                }
            }
        }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::BorderAdd( LE * a_field, TTypesOfLandscapeElement a_type ) {
    int x = a_field->GetValidX();
    int y = a_field->GetValidY();
  if ( ( x == -1 ) || ( y == -1 ) ) {
    g_msg->Warn( WARN_BUG, "Landscape::BorderAdd(): Uninitialized border coordinate!", "" );
    exit( 1 );
  }
  LE * border = NewElement(a_type);
  a_field->SetBorder( border );
  m_polymapping[ hb_first_free_poly_num ] = (int) m_elems.size();
  m_elems.resize( m_elems.size() + 1 );
  m_elems[ m_elems.size() - 1 ] = border;
  border->SetPoly( hb_first_free_poly_num++ );
  border->SetArea( 0.0 );
  BorderScan(a_field, g_map_le_borderwidth.value());
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::BorderScan( LE * a_field, int a_width ) 
{
    LE * border = a_field->GetBorder(); // border is the a border object
    int fieldpoly = a_field->GetPoly(); // fieldpoly is the polygon number
    int borderpoly = border->GetPoly(); // borderpoly is the polygon number
    int borderindex = m_polymapping[ borderpoly ]; // borderindex is the elems index for the border
    int fieldindex = m_polymapping[ fieldpoly ]; // fieldindex is the elems index
    int test = m_land->Get(a_field->GetCentroidX(), a_field->GetCentroidY());
    if (test != fieldindex)
    {
        g_msg->Warn("Landscape::BorderScan - Border Scan centroid does not return correct polygon index. Index :", fieldindex); 
        g_msg->Warn(" Returned ", test);
        exit(0);
    }
    int notforever = 50000;
    vector<APoint> listoflocs;
    for (int wid = 0; wid < a_width; wid++)
    {
        notforever = 50000;
        // These two will be modified through pointer operations in BorderStep().
        APoint coord(a_field->GetCentroidX(), a_field->GetCentroidY());
        // Find the first edge cell
        AxisLoop(fieldindex, &coord, random(8));
        while (--notforever > 0)
        {
            // Check if this position should be made into a border.
            if (BorderTest(fieldindex, -99, coord.m_x, coord.m_y))
            {
                // Add this pixel to the border element in the big map, but using a code for later replacement.
                m_land->Put(coord.m_x, coord.m_y, -99); // this puts the elems index into our map in memory
                listoflocs.push_back(coord);
                a_field->AddArea(-1.0);
                if (l_map_exit_on_zero_area.value() && (a_field->GetArea()<1))
                {
                    char polynum[20];
                    sprintf(polynum, "%d", a_field->GetPoly()); g_msg->Warn(WARN_FILE, "Landscape::BorderScan(): Polygon reached zero area " "when adding border. Poly num: ", polynum);
                    exit(1);
                }
                border->AddArea(1.0);
                border->SetMapValid(true);
            }
            // Step to next coordinate. Quit when done.
            if (!BorderStep(fieldindex, -99, &coord))
            {
                break;
            }
        }
        for (std::vector<APoint>::iterator it = listoflocs.begin(); it != listoflocs.end(); ++it)
        {
            m_land->Put((*it).m_x, (*it).m_y, borderindex);
        }
        listoflocs.clear();
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::StepOneValid( int a_polyindex, int a_x, int a_y, int a_step )
{
  int index;
  int x_add[ 8 ] = { 1*a_step, 1*a_step, 0, -1*a_step, -1*a_step, -1*a_step, 0, 1*a_step };
  int y_add[ 8 ] = { 0, -1*a_step, -1*a_step, -1*a_step, 0, 1*a_step, 1*a_step, 1*a_step };
  int width = m_land->MapWidth();
  int height = m_land->MapHeight();
  // Scan anti-clockwise from center pixel coordinate.
  for ( unsigned int i = 0; i < 8; i++ ) {
    if ( ( a_x + x_add[ i ] < width ) && ( a_x + x_add[ i ] >= 0 ) && ( a_y + y_add[ i ] < height ) && ( a_y + y_add[ i ] >= 0 ) )
    {
        index = m_land->Get( a_x + x_add[ i ], a_y + y_add[ i ] );
        if (  index == a_polyindex )
        {
            m_elems[a_polyindex]->SetValidXY(a_x + x_add[ i ], a_y + y_add[ i ]);
            return true;
        }
    }
  }
  return false;
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::BorderTest( int a_fieldindex, int a_borderindex, int a_x, int a_y )
{
  int index;
  int x_add[ 8 ] = { 1, 1, 0, -1, -1, -1, 0, 1 };
  int y_add[ 8 ] = { 0, -1, -1, -1, 0, 1, 1, 1 };
  int width = m_land->MapWidth();
  int height = m_land->MapHeight();
  // Scan anti-clockwise from center pixel coordinate.
  for ( unsigned int i = 0; i < 8; i++ ) {
    if ( ( a_x + x_add[ i ] >= width ) || ( a_x + x_add[ i ] < 0 ) || ( a_y + y_add[ i ] >= height )
         || ( a_y + y_add[ i ] < 0 ) ) {
           return true;
    }
    //continue;
    index = m_land->Get( a_x + x_add[ i ], a_y + y_add[ i ] );
    if ( ( index != a_fieldindex ) && ( index != a_borderindex ) )
    {
        return true;
        // Test removed 1/07/2014 CJT
      //if ( BorderNeed( m_elems[ index ]->GetElementType() ) ) return true;
      //else return false;
    }
  }
  return false;
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::BorderStep(int a_fieldindex, int /* a_borderindex */, int * a_x, int * a_y) {
    int index;
    int x_add[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
    int y_add[8] = { 0, -1, -1, -1, 0, 1, 1, 1 };
    int width = m_land->MapWidth();
    int height = m_land->MapHeight();
    int i = 7, counter = 8;
    bool running = true;
    // First scan for another pixel that belongs to this field.
    while (running)
    {
        if (!((*a_x) + x_add[i] >= width) && !((*a_x) + x_add[i] < 0) && !((*a_y) + y_add[i] >= height) && !((*a_y) + y_add[i] < 0))
        {
            index = m_land->Get((*a_x) + x_add[i], (*a_y) + y_add[i]);
            if (index == a_fieldindex)
            {
                // Found the first field pixel while scanning around always
                // in the same direction.
                running = false;
            }
        }
        if (--i < 0) {
            // Didn't find any of our pixels. We are in a blind alley. Exit
            // gracefully.
            return false; // Signal done scanning this field.
        }
    }

    // Now scan around from our present facing direction and find the border
    // (if any).
    while (--counter)
    {
        if (!((*a_x) + x_add[i] >= width) && !((*a_x) + x_add[i] < 0) && !((*a_y) + y_add[i] >= height) && !((*a_y) + y_add[i] < 0))
        {
            index = m_land->Get((*a_x) + x_add[i], (*a_y) + y_add[i]);
            if (index == a_fieldindex)
            {
                if (--i < 0) i = 7;
                continue;
            }
        }

        // Aha! This pixel is not ours. Step one step in the
        // opposite(!) direction. If that pixel is ours, then
        // modify hotspot coordinates and exit.
        if (++i > 7) i = 0;
        if (!((*a_x) + x_add[i] + 1 > width) && !((*a_x) + x_add[i] < 0) && !((*a_y) + y_add[i] + 1 > height) &&
            !((*a_y) + y_add[i] < 0) && (m_land->Get((*a_x) + x_add[i], (*a_y) + y_add[i]) == a_fieldindex))
        {
            (*a_x) += x_add[i];
            (*a_y) += y_add[i];
            return true;
        }
    }
    return false;
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::BorderStep(int a_fieldindex, int /* a_borderindex */, APoint* a_coord) {
    int index;
    int x_add[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
    int y_add[8] = { 0, -1, -1, -1, 0, 1, 1, 1 };
    int width = m_land->MapWidth();
    int height = m_land->MapHeight();
    int i = 7, counter = 8;
    bool running = true;
    // First scan for another pixel that belongs to this field.
    while (running)
    {
        if (!((a_coord->m_x) + x_add[i] >= width) && !((a_coord->m_x) + x_add[i] < 0) && !((a_coord->m_y) + y_add[i] >= height) && !((a_coord->m_y) + y_add[i] < 0))
        {
            index = m_land->Get((a_coord->m_x) + x_add[i], (a_coord->m_y) + y_add[i]);
            if (index == a_fieldindex)
            {
                // Found the first field pixel while scanning around always
                // in the same direction.
                running = false;
            }
        }
        if (--i < 0) {
            // Didn't find any of our pixels. We are in a blind alley. Exit
            // gracefully.
            return false; // Signal done scanning this field.
        }
    }

    // Now scan around from our present facing direction and find the border
    // (if any).
    while (--counter)
    {
        if (!((a_coord->m_x) + x_add[i] >= width) && !((a_coord->m_x) + x_add[i] < 0) && !((a_coord->m_y) + y_add[i] >= height) && !((a_coord->m_y) + y_add[i] < 0))
        {
            index = m_land->Get((a_coord->m_x) + x_add[i], (a_coord->m_y) + y_add[i]);
            if (index == a_fieldindex)
            {
                if (--i < 0) i = 7;
                continue;
            }
        }

        // Aha! This pixel is not ours. Step one step in the
        // opposite(!) direction. If that pixel is ours, then
        // modify hotspot coordinates and exit.
        if (++i > 7) i = 0;
        if (!((a_coord->m_x) + x_add[i] + 1 > width) && !((a_coord->m_x) + x_add[i] < 0) && !((a_coord->m_y) + y_add[i] + 1 > height) &&
            !((a_coord->m_y) + y_add[i] < 0) && (m_land->Get((a_coord->m_x) + x_add[i], (a_coord->m_y) + y_add[i]) == a_fieldindex))
        {
            (a_coord->m_x) += x_add[i];
            (a_coord->m_y) += y_add[i];
            return true;
        }
    }
    return false;
}
//-----------------------------------------------------------------------------------------------------------------------------
void Landscape::CreatePondList()
{
    for (unsigned int i = 0; i < m_elems.size(); i++)
    {
        if (m_elems[i]->GetElementType() == tole_Pond) {
            m_PondIndexList.push_back(i);
            m_PondRefsList.push_back(m_elems[i]->GetPoly());
        }
    }
}

int Landscape::SupplyRandomPondIndex() { 
    return m_PondIndexList[int(g_rand_uni()*m_PondIndexList.size())]; 
}

int Landscape::SupplyRandomPondRef() { 
    if (m_PondIndexList.size()>0) return m_PondRefsList[int(g_rand_uni()*m_PondIndexList.size())]; 
    return -1;
}


// Unsprayed Margin Code....
void Landscape::UnsprayedMarginAdd( LE * a_field ) {
  int x = a_field->GetValidX();
  int y = a_field->GetValidY();
  if ( ( x == -1 ) || ( y == -1 ) ) {
    // Tripping this probably means it is not a field
    g_msg->Warn( WARN_BUG, "Landscape::UnsprayedMarginAdd(): Uninitialized border coordinate!", "" );
    exit( 1 );
  }
  LE * umargin = NewElement( tole_UnsprayedFieldMargin );
  m_polymapping[ hb_first_free_poly_num ] =  (int) m_elems.size();
  m_elems.resize( m_elems.size() + 1 );
  m_elems[ m_elems.size() - 1 ] = umargin;
  a_field->SetUnsprayedMarginPolyRef( hb_first_free_poly_num );
  umargin->SetPoly( hb_first_free_poly_num++ );
  umargin->SetArea( 0.0 );

  for ( int q = 0; q < l_map_umargin_width.value(); q++ )
    UnsprayedMarginScan( a_field, q + 1 );
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::UnsprayedMarginScan( LE * a_field, int a_width ) {
  LE * umargin = g_landscape_p->SupplyLEPointer( a_field->GetUnsprayedMarginPolyRef() );
  int fieldpoly = a_field->GetPoly();
  int borderpoly = umargin->GetPoly();
  int borderindex = m_polymapping[ borderpoly ];
  int fieldindex = m_polymapping[ fieldpoly ];
  int notforever = 5000;

  // These two will be modified through pointer operations
  // in BorderStep().
  int x = a_field->GetValidX();
  int y = a_field->GetValidY();
  // Now the problem is that GetValid does not always return a valid co-ord!
  // so we need to search for one
  if ( !FindValidXY( fieldindex, x, y ) ) return;

  while ( --notforever ) {
    // Check if this position should be made into a border.
    if ( UMarginTest( fieldindex, borderindex, x, y, a_width ) ) {
      // Add this pixel to the border element in the big map.
      m_land->Put( x, y, borderindex );
      a_field->AddArea( -1.0 );
      umargin->AddArea( 1.0 );
    };
    // Step to next coordinate. Quit when done.
    if ( !BorderStep( fieldindex, borderindex, & x, & y ) )
      return;
  }
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::UMarginTest( int a_fieldindex, int a_marginindex, int a_x, int a_y, int a_width ) {
  int index;
  int x_add[ 8 ] = { 1*a_width, 1*a_width, 0, -1*a_width, -1*a_width, -1*a_width, 0, 1*a_width };
  int y_add[ 8 ] = { 0, -1*a_width, -1*a_width, -1*a_width, 0, 1*a_width, 1*a_width, 1*a_width };
  int width = m_land->MapWidth();
  int height = m_land->MapHeight();
  // Scan anti-clockwise from center pixel coordinate.
  for ( unsigned int i = 0; i < 8; i++ ) {
    if ( ( a_x + x_add[ i ] >= width ) || ( a_x + x_add[ i ] < 0 ) || ( a_y + y_add[ i ] >= height )
         || ( a_y + y_add[ i ] < 0 ) ) {
           return true;
    }
    //continue;
    index = m_land->Get( a_x + x_add[ i ], a_y + y_add[ i ] );
    if ( ( index != a_fieldindex ) && ( index != a_marginindex ) ) return true;
  }
  return false;
}

bool Landscape::FindValidXY( int a_field, int & a_x, int & a_y ) {
  // From a hopefully sensible starting point this method scans in the
  // 8 directions to find a good valid x and y matching a_field
  int x_add[ 8 ] = { 1, 1, 0, -1, -1, -1, 0, 1 };
  int y_add[ 8 ] = { 0, -1, -1, -1, 0, 1, 1, 1 };
  int index;
  int nx, ny;
  int width = m_land->MapWidth();
  int height = m_land->MapHeight();
  // Assume it has to within 100m
  for ( int i = 0; i < 100; i++ ) {
    for ( int l = 0; l < 8; l++ ) {
      nx = a_x + x_add[ l ] * i;
      ny = a_y + y_add[ l ] * i;
      if ( ( nx < width ) && ( nx >= 0 ) && ( ny < height ) && ( ny >= 0 ) ) {
        index = m_land->Get( nx, ny );
        if ( index == a_field ) {
          a_x = a_x + x_add[ l ] * i;
          a_y = a_y + y_add[ l ] * i;
          return true;
        }
      }
    }
  }
  return false;
}
//-----------------------------------------------------------------------------------------------------------------------------
// ...to here..S

// **************************************************************************************
// **************************************************************************************

void Landscape::AddBeetleBanks( TTypesOfLandscapeElement a_tole ) {
    int BBs=0;
    int tx1 = cfg_BeetleBankMinX.value();
    int tx2 = cfg_BeetleBankMaxX.value();
    int ty1 = cfg_BeetleBankMinY.value();
    int ty2 = cfg_BeetleBankMaxY.value();
    bool doit = false;
    unsigned sz=(unsigned) m_elems.size();
    for (unsigned i=0; i<sz; i++)
    {
        if (m_elems[ i ]->GetElementType() == tole_Field)
        {
            doit = false;
            int cx = m_elems[ i ]->GetCentroidX();
            int cy = m_elems[ i ]->GetCentroidY();
            if (!cfg_BeetleBankInvert.value())
            {
                if ((cx >= tx1) && (cy >= ty1) && (cx <= tx2) && (cy <= ty2))
                {
                    doit = true;
                }
            }
            else if ((cx < tx1) || (cy < ty1) || (cx > tx2) || (cy > ty2))
            {
                doit = true;
            }
            if (doit)
            {
                if (random(100)<cfg_BeetleBankChance.value())
                {
                    if (BeetleBankPossible( m_elems[ i ], a_tole) ) BBs++;
                }
            }
        }
    }
    char str[25];
    sprintf(str,"%d",BBs);
    g_msg->Warn( WARN_MSG, "Landscape::AddBeetleBanks(): BeetleBanks successfully added:", str );
}
//-----------------------------------------------------------------------------------------------------------------------------
// **************************************************************************************

bool Landscape::BeetleBankPossible( LE * a_field, TTypesOfLandscapeElement a_tole ) {
    int farea=(int)a_field->GetArea();
    if (farea<10000) return false;
    int cx=a_field->GetCentroidX();
    int cy=a_field->GetCentroidY();
    // The centroid is the only estimate we have (and it at least should be in the field).
    // So start here and find the centre
    if (!FindFieldCenter(a_field, &cx, &cy)) return false;
    // now get the alignment
    int length=0;
    int alignment=FindLongestAxis(&cx, &cy, &length);
    // reduce length by 20%
    length=int(length*0.8);
    int area=2*length*cfg_BeetleBankWidth.value(); // 12m wide fixed size
    if (area>(farea*cfg_BeetleBankMaxArea.value())) return false;
    // Must be small engough so lets draw it
    BeetleBankAdd(cx, cy, alignment, length , a_field, a_tole);
    return true;
}
//-----------------------------------------------------------------------------------------------------------------------------

bool Landscape::FindFieldCenter(LE* a_field, int* x, int* y) {
    // Start at x,y
    // works by selecting the point that is a mean of the co-ords of the centers of 4 axes from this point that are in the field.
    // Then do it again, and again until we don't move more than 1m or we have tried too many times
    int ourpoly=SupplyPolyRef(*(x),*(y));
    if (ourpoly!=a_field->GetPoly()) return false;
    int centers[2][8];
    int tries=0;
    int diff=999;
    int x1=*(x);
    int y1=*(y);
    int centreX=x1;
    int centreY=y1;
    // NB we might escape without bounds checking here because the polygon number does not wrap round - will only ever be a problem if we go SimX+1,SimY+1
    while ((diff>1) & (tries++<100)) {
        for (unsigned v=0; v<4; v++) {
            x1=centreX;
            y1=centreY;
            AxisLoop(ourpoly, &x1, &y1, v);
            centers[0][v]=x1-m_x_add[v];
            centers[1][v]=y1-m_y_add[v];
            x1=centreX;
            y1=centreY;
            AxisLoop(ourpoly, &x1, &y1, v+4);
            centers[0][v+4]=x1-m_x_add[v+4];
            centers[1][v+4]=y1-m_y_add[v+4];
//          centreX+=((centers[0][v]+x1-m_x_add[v+4])/2);
//          centreY+=((centers[1][v]+y1-m_y_add[v+4])/2);
        }
        int oldx=centreX;
        int oldy=centreY;
        centreX=0;
        centreY=0;
        for (int h=0; h<8; h++) {
            centreX+=centers[0][h];
            centreY+=centers[1][h];
        }
        centreX/=8;
        centreY/=8;
        diff=abs(oldx-centreX)+abs(oldy-centreY);
    }
    *(x)=centreX;
    *(y)=centreY;
    int tourpoly=SupplyPolyRef(*(x),*(y));
    if (tourpoly!=ourpoly) {
        return false; // can happen eg if there is a pond in the middle of the field
    }

    return true;
}
//-----------------------------------------------------------------------------------------------------------------------------

int Landscape::FindLongestAxis(int* a_x, int* a_y, int* a_length) 
{
    int ourpoly=SupplyPolyRef(*(a_x),*(a_y));
    int dist[4];
    int distx[8];
    int disty[8];
    int found = -1;
    *(a_length) = 0;
    int dx[8];
    int dy[8];
    int fx[8];
    int fy[8];
    for (unsigned v=0; v<8; v++) 
    {
        int x1=*(a_x);
        int y1=*(a_y);
        AxisLoop(ourpoly, &x1, &y1, v);
        x1 -= m_x_add[v];
        y1 -= m_y_add[v];
        dx[v] = abs(*(a_x)-x1);
        dy[v] = abs(*(a_y)-y1);
        fx[v] = x1;
        fy[v] = y1;
        distx[v] = dx[v];
        disty[v] = dy[v];
    }
    for (int di = 0; di < 4; di++)
    {
        int dx = distx[di] + distx[di + 4];
        int dy = disty[di] + disty[di + 4];
        if (dx == 0) dist[di] = dy; else dist[di] = dx;
        if (dist[di] > *(a_length))
        {
            found = di;
            *(a_length) = dist[di];
        }
    }
    if (found == -1) return 0;
    // Now need to find the middle of the axis.
    int l = (*(a_length) / 2);
    if (fx[found] > fx[found + 4]) *(a_x) = fx[found + 4] + m_x_add[found] * l;  else *(a_x) = fx[found + 4] - m_x_add[found + 4] * l;
    if (fy[found] > fy[found + 4]) *(a_y) = fy[found + 4] + m_y_add[found] * l;  else *(a_y) = fy[found + 4] - m_y_add[found + 4] * l;

    return found;
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::AxisLoop(int a_polyindex, APoint* a_cor, int a_axis) {
    int ap1 = a_polyindex;
    while (ap1 == a_polyindex)
    {
        a_cor->m_x += m_x_add[a_axis];
        a_cor->m_y += m_y_add[a_axis];
        if (a_cor->m_x >= m_width - 1) { a_cor->m_x = m_width - 1; return; }
        if (a_cor->m_y >= m_height - 1) { a_cor->m_y = m_height - 1; return; }
        if (a_cor->m_x <= 0) { a_cor->m_x = 0; return; }
        if (a_cor->m_y <= 0) { a_cor->m_y = 0; return; }
        ap1 = m_land->Get(a_cor->m_x, a_cor->m_y); // NB this returns the m_elemens index not the polyref (ChangeMapMapping has been called by here)
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::AxisLoopLtd(int a_polyindex, APoint* a_cor, int a_axis, int a_limit) {
    int ap1 = a_polyindex;
    int count = 0;
    while (ap1 == a_polyindex && count<a_limit)
    {
        a_cor->m_x += m_x_add[a_axis];
        a_cor->m_y += m_y_add[a_axis];
        if (a_cor->m_x >= m_width - 1) { a_cor->m_x = m_width - 1; return; }
        if (a_cor->m_y >= m_height - 1) { a_cor->m_y = m_height - 1; return; }
        if (a_cor->m_x <= 0) { a_cor->m_x = 0; return; }
        if (a_cor->m_y <= 0) { a_cor->m_y = 0; return; }
        ap1 = m_land->Get(a_cor->m_x, a_cor->m_y); // NB this returns the m_elemens index not the polyref (ChangeMapMapping has been called by here)
        count++;
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::AxisLoop(int a_polyindex, int* a_x, int* a_y, int a_axis) {
    int ap1 = a_polyindex;
    while (ap1 == a_polyindex)
    {
        *(a_x) += m_x_add[a_axis];
        *(a_y) += m_y_add[a_axis];
        // Before we try to get a polyindex from the map, check we are still on the world
        if (*(a_x) < 0) 
        {
            return;
        }
        if (*(a_y) < 0)
        {
            return;
        }
        if (*(a_x) >= m_width) 
        {
            return;
        }
        if (*(a_y) >= m_height) 
        {
            return;
        }
        // OK still in the map, get the polyindex
        ap1 = m_land->Get((*a_x), (*a_y)); // NB this returns the m_elemens index not the polyref (ChangeMapMapping has been called by here)
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::BeetleBankAdd(int a_x, int a_y, int a_angle, int a_length, LE* a_field, TTypesOfLandscapeElement a_tole ) {
    // Need to get a new number
    ++m_LargestPolyNumUsed;
    // Make the new landscape element
    LE * BeetleBank;
    switch (a_tole)
    {
    case tole_MownGrass:
        BeetleBank = NewElement( tole_MownGrass );
        BeetleBank->SetALMaSSEleType( g_letype->BackTranslateEleTypes( tole_MownGrass ) );
        break;
    case tole_PermanentSetaside:
        BeetleBank = NewElement( tole_PermanentSetaside );
        BeetleBank->SetALMaSSEleType( g_letype->BackTranslateEleTypes( tole_PermanentSetaside ) );
        break;
    case tole_BeetleBank:
    default:
        BeetleBank = NewElement( tole_BeetleBank );
        BeetleBank->SetALMaSSEleType( g_letype->BackTranslateEleTypes( tole_BeetleBank ) );
    }
    BeetleBank->SetVegPatchy(true);
    m_polymapping[ m_LargestPolyNumUsed ] =  (int) m_elems.size();
    m_elems.resize( m_elems.size() + 1 );
    m_elems[ m_elems.size() - 1 ] = BeetleBank;
    BeetleBank->SetPoly( m_LargestPolyNumUsed );
    // write lengthx12m to the map at alignment angle
    int area=0;
    int angle2=0;
    int width=cfg_BeetleBankWidth.value();
    if (a_angle==0) angle2=2;
    int start=(int)(a_length*0.1);
    for (int i=start; i<a_length; i++) {
        for (int w=0-width; w<width; w++) {
            int tx=w*m_x_add[angle2];
            int ty=w*m_y_add[angle2];
            m_land->Put( tx+a_x+i*m_x_add[a_angle], ty+a_y+i*m_y_add[a_angle], (int) m_elems.size() - 1 );
            m_land->Put( tx+a_x-i*m_x_add[a_angle], ty+a_y-i*m_y_add[a_angle], (int) m_elems.size() - 1 );
            area+=2;
            a_field->AddArea( -2.0 );

        }
    }
    BeetleBank->SetArea( double(area) );
    BeetleBank->SetValidXY( a_x+start*m_x_add[a_angle], a_y+start*m_y_add[a_angle] );
    BeetleBank->SetMapValid(true);

}
//-----------------------------------------------------------------------------------------------------------------------------

// **************************************************************************************
/* End beetle bank addition code */
// **************************************************************************************/


void Landscape::SkylarkEvaluation(SkTerritories* a_skt) {
    for (unsigned i=0; i<m_elems.size(); i++) {
        a_skt->PreCachePoly(m_elems[i]->GetPoly());
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::RodenticidePredatorsEvaluation(RodenticidePredators_Population_Manager* a_rppm)
{
    for (unsigned i=0; i<m_elems.size(); i++) {
        a_rppm->PreCachePoly(m_elems[i]->GetPoly());
    }
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::CalculateCentroids()
{
    cout << "In Centroid Calculations" << endl;
    // For each polygon
    for (int p = 0; p< (int)m_elems.size(); p++)
    {
        // Calcuate the actual centre
        int x1 = m_elems[p]->GetMinX();
        int y1 = m_elems[p]->GetMinY();
        int x2 = m_elems[p]->GetMaxX();
        int y2 = m_elems[p]->GetMaxY();
        int midx = (x1 + x2) / 2;
        int midy = (y1 + y2) / 2;
        // Now from midx & midy we move outwards in concentric circles until we find a location that matches our polyref.
        int polyindex = p; // Change mapmapping has been called by now, so the map contains m_elems indices.
        CentroidSpiralOut(polyindex, midx, midy);
        // Now we want to be sure that we are in the middle of the polygon not on the edge. This is tricky for complex shaped polygons, 
        // but we have a stab at it by using the FindLongestAxis method. This puts us in the centre of the longest axis in 8 directions
        // from this point
        int l;
        FindLongestAxis(&midx, &midy, &l);
        m_elems[p]->SetCentroid(midx, midy);
    }
    BuildingDesignationCalc();
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::CentroidSpiralOut(int a_polyref, int &a_x, int &a_y)
{
    if (SupplyPolyRefIndex(a_x, a_y) == a_polyref) return; // Found it so return
    // Otherwise its not found so we need to start to spiral out
    int loop = 1;
    int sx = a_x;
    int sy = a_y;
    do {
        a_y = sy - loop;
        for (int i = 0 - loop; i <= loop; i++)
        {
            a_x = sx + i;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
                return; // Found it so return
        }
        a_y = sy + loop;
        for (int i = 0 - loop; i <= loop; i++)
        {
            a_x = sx + i;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
                return; // Found it so return
        }
        a_x = sx + loop;
        for (int j = 0 - (loop - 1); j< loop; j++)
        {
            a_y = sy + j;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
                return; // Found it so return
        }
        a_x = sx - loop;
        for (int j = 0 - (loop - 1); j< loop; j++)
        {
            a_y = sy + j;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
                return; // Found it so return
        }
        loop++;
    } while (loop<m_width); // This stopping rule should hopefully not be needed, it is set very high.
    g_msg->Warn("Landscape::CentroidSpiralOut: Failure of centroid main loop. Looking for polygon index ",a_polyref);
    a_x = m_elems[a_polyref]->GetMinX();
    a_y = m_elems[a_polyref]->GetMinY();
}
//-----------------------------------------------------------------------------------------------------------------------------
void Landscape::CentroidSpiralOutBlocks(int a_polyref, int &a_x, int &a_y)
{
    if (SupplyPolyRefIndex(a_x, a_y) == a_polyref) return; // Found it so return
    // Otherwise its not found so we need to start to spiral out until we find a 10x10m block of our field.
    int loop = 1;
    int sx = a_x;
    int sy = a_y;
    do {
        a_y = sy - loop;
        for (int i = 0 - loop; i <= loop; i++)
        {
            a_x = sx + i;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
            {
                APoint pt(a_x, a_y);
                double dists = -1;
                double best = -1;
                int ind = -1;
                for (int d = 0; d < 8; d++)
                {
                    AxisLoopLtd(a_polyref, &pt, d,10);
#ifdef __BORLANDC__
                    dists = sqrt(double((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y)));
#else
                    //dists = sqrt((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y));
                    dists = (abs(a_x - pt.m_x) + abs(a_y - pt.m_y));
#endif
                    if (dists > best) { best = dists; ind = d; }
                }
                if (ind > -1)
                {
                    //                  best /= 2;
                    if (best >= 10)
                    {
                        a_x += m_x_add[ind] * 10;
                        a_y += m_y_add[ind] * 10;
                        return;
                    }
                }
                return; // Found it so return even without the 10 distance
            }
        }
        a_y = sy + loop;
        for (int i = 0 - loop; i <= loop; i++)
        {
            a_x = sx + i;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
            {
                APoint pt(a_x, a_y);
                double dists = -1;
                double best = -1;
                int ind = -1;
                for (int d = 0; d < 8; d++)
                {
                    AxisLoopLtd(a_polyref, &pt, d,10);
#ifdef __BORLANDC__
                    dists = sqrt(double((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y)));
#else
                    //dists = sqrt((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y));
                    dists = (abs(a_x - pt.m_x) + abs(a_y - pt.m_y));
#endif
                    if (dists > best) { best = dists; ind = d; }
                }
                if (ind > -1)
                {
                    //                  best /= 2;
                    if (best >= 10)
                    {
                        a_x += m_x_add[ind] * 10;
                        a_y += m_y_add[ind] * 10;
                        return;
                    }
                }
                return; // Found it so return even without the 10 distance
            }
        }
        a_x = sx + loop;
        for (int j = 0 - (loop - 1); j< loop; j++)
        {
            a_y = sy + j;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
            {
                APoint pt(a_x, a_y);
                double dists = -1;
                double best = -1;
                int ind = -1;
                for (int d = 0; d < 8; d++)
                {
                    AxisLoopLtd(a_polyref, &pt, d,10);
#ifdef __BORLANDC__
                    dists = sqrt(double((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y)));
#else
                    //dists = sqrt((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y));
                    dists = (abs(a_x - pt.m_x) + abs(a_y - pt.m_y));
#endif
                    if (dists > best) { best = dists; ind = d; }
                }
                if (ind > -1)
                {
                    //                  best /= 2;
                    if (best >= 10)
                    {
                        a_x += m_x_add[ind] * 10;
                        a_y += m_y_add[ind] * 10;
                        return;
                    }
                }
                return; // Found it so return even without the 10 distance
            }
        }
        a_x = sx - loop;
        for (int j = 0 - (loop - 1); j< loop; j++)
        {
            a_y = sy + j;
            CorrectCoords(a_x, a_y);
            if (SupplyPolyRefIndex(a_x, a_y) == a_polyref)
            {
                APoint pt(a_x, a_y);
                double dists = -1;
                double best = -1;
                int ind = -1;
                for (int d = 0; d < 8; d++)
                {
                    AxisLoopLtd(a_polyref, &pt, d,10);
#ifdef __BORLANDC__
                    dists = sqrt(double((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y)));
#else
                    //dists = sqrt((a_x - pt.m_x)*(a_x - pt.m_x) + (a_y - pt.m_y) * (a_y - pt.m_y));
                    dists = (abs(a_x - pt.m_x) + abs(a_y - pt.m_y));
#endif
                    if (dists > best) { best = dists; ind = d; }
                }
                if (ind > -1)
                {
                    //                  best /= 2;
                    if (best >= 10)
                    {
                        a_x += m_x_add[ind] * 10;
                        a_y += m_y_add[ind] * 10;
                        return;
                    }
                }
                return; // Found it so return even without the 10 distance
            }
        }
        loop++;
    } while (loop<m_width); // This stopping rule should hopefully not be needed, it is set very high.
    exit(0);
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::DumpCentroids()
{
    ofstream centroidfile("PolygonCentroids.txt", ios::out);
    centroidfile<<"Polyref"<<'\t'<<"CX"<<'\t'<<"CY"<<'\t'<<"Type"<<'\t'<<"Area"<<'\t'<<"Country Designation"<<endl;
    for (int p = 0; p< (int)m_elems.size(); p++)
    {
        centroidfile<<m_elems[p]->GetPoly()<<'\t'<<m_elems[p]->GetCentroidX()<<'\t'<<m_elems[p]->GetCentroidY()<<'\t'<<m_elems[p]->GetElementType()<<'\t'<<m_elems[p]->GetArea()<<'\t'<<m_elems[p]->GetCountryDesignation()<<endl;
    }
    centroidfile.close();
}
//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::BuildingDesignationCalc()
{
    for (int p = 0; p< (int)m_elems.size(); p++)
    {
        TTypesOfLandscapeElement tole = m_elems[p]->GetElementType();
        if ( tole == tole_Building)
        {
            int cx = m_elems[p]->GetCentroidX();
            int cy = m_elems[p]->GetCentroidY();
            int near = 0;
            for (int j = 0; j< (int)m_elems.size(); j++)
            {
                if (m_elems[j]->GetElementType() == tole_Building)
                {
                    int nx = m_elems[j]->GetCentroidX();
                    int ny = m_elems[j]->GetCentroidY();
                    int dx =abs(cx-nx);
                    int dy =abs(cy-ny);
                    if ((dx < cfg_mintownbuildingdistance.value()) && (dy < cfg_mintownbuildingdistance.value())) near++;
                    if (near > cfg_mintownbuildingdistance.value()) break;
                }
            }
            if (near <= cfg_mintownbuildingnumber.value()) m_elems[p]->SetCountryDesignation(1); // Not enough buildings close by, so it is a country building
            else m_elems[p]->SetCountryDesignation(0);
        }
        else if (tole == tole_YoungForest) 
        {
            m_elems[p]->SetCountryDesignation(2);
        }
        else if ((tole == tole_DeciduousForest) || ( tole == tole_MixedForest) || ( tole == tole_ConiferousForest ) ) m_elems[p]->SetCountryDesignation(3);
            
    }
}
//-----------------------------------------------------------------------------------------------------------------------------


void Landscape::RecordGooseNumbers(int a_polyref, int a_number)
{
    int day = SupplyDayInYear();
    m_elems[m_polymapping[a_polyref]]->SetGooseNos(a_number, day);
}

//-----------------------------------------------------------------------------------------------------------------------------

void Landscape::RecordGooseNumbersTimed(int a_polyref, int a_number)
{
    int day = SupplyDayInYear();
    m_elems[m_polymapping[a_polyref]]->SetGooseNosTimed(a_number, day);
}

//-----------------------------------------------------------------------------------------------------------------------------
void Landscape::RecordGooseSpNumbers(int a_polyref, int a_number, GooseSpecies a_goose) {
      int day = SupplyDayInYear();
      m_elems[m_polymapping[a_polyref]]->SetGooseSpNos(a_number, day, a_goose);
  }
  //-----------------------------------------------------------------------------------------------------------------------------
  void Landscape::RecordGooseSpNumbersTimed(int a_polyref, int a_number, GooseSpecies a_goose) {
      int day = SupplyDayInYear();
      m_elems[m_polymapping[a_polyref]]->SetGooseSpNosTimed(a_number, day, a_goose);
  }
  //-----------------------------------------------------------------------------------------------------------------------------
  int Landscape::GetGooseNumbers(int a_polyref) {
      return m_elems[m_polymapping[a_polyref]]->GetGooseNos();
  }
  //-----------------------------------------------------------------------------------------------------------------------------

  int Landscape::GetGooseNumbers(int a_x, int a_y)
  {
      return m_elems[m_land->Get(a_x, a_y)]->GetGooseNos();
  }
  //-----------------------------------------------------------------------------------------------------------------------------

double Landscape::GetHareFoodQuality(int a_polygon)  // actually only works for a single square, but who cares, a polygon is uniform
{
          double digest;
          TTypesOfLandscapeElement habitat = SupplyElementType(a_polygon);
          switch (habitat) {
              // Impossible stuff
          case tole_Building:
          case tole_Pond:
          case tole_Freshwater:
          case tole_River:
          case tole_Saltwater:
          case tole_Coast:
          case tole_BareRock:
          case tole_ConiferousForest:
          case tole_DeciduousForest:
          case tole_MixedForest:
          case tole_SmallRoad:
          case tole_LargeRoad:
          case tole_ActivePit:
          case tole_UrbanNoVeg:
          case tole_UrbanPark:
          case tole_SandDune:
          case tole_Copse:
          case tole_Stream:
          case tole_MetalledPath:
          case tole_Carpark:
          case tole_FishFarm:
              //            EnergyBalance(activity_Foraging, 100); // This is a bug - it penalises for foraging in impossible areas - not intended but not found until after parameter fitting! Removed 28/07/2014
              return 0.0;

              // Questionable stuff
          case tole_RiversidePlants:
          case tole_RiversideTrees:
          case tole_Garden:
          case tole_Track:
          case tole_StoneWall:
          case tole_Hedges:
          case tole_Marsh:
          case tole_PitDisused:
          case tole_RoadsideVerge:
          case tole_Railway:
          case tole_Scrub:
          case tole_AmenityGrass:
          case tole_Parkland:
          case tole_BuiltUpWithParkland:
          case tole_Churchyard:
          case tole_HeritageSite:
              return 0.25; // was 0.25  being half of access to low digestability stuff
              //        case tole_MownGrass:
              //            digest = 0.8;  // Added 28/07/2014 this is a way to compensate for the lack of choice when foraging, i.e. the whole area  is assumed to be foraged equally.
          case tole_Wasteland:
          case tole_IndividualTree:
          case tole_WoodyEnergyCrop:
          case tole_PlantNursery:
          case tole_Pylon:
          case tole_WindTurbine:
          case tole_WoodlandMargin:
          case tole_Vildtager:
          default:
              digest = SupplyVegDigestability(a_polygon);
          }
#ifdef __Perfectfood
          return 0.8;
#else
#ifdef __YEARLYVARIABLEFOODQUALITY
          digest *= m_OurPopulationManager->m_GoodYearBadYear;
#endif
          double veg_height;
          double access = 1.0;
          // double grazedreduction[4] = { 1.0, 0.75, 0.5, 0.25 };
          double grazedreduction[4] = { 1.0, 0.8, 0.2, 0.05 };
          veg_height = SupplyVegHeight(a_polygon);
          double weeds = SupplyWeedBiomass(a_polygon);
          if ((veg_height <= 0) && (weeds < 0.1)) return 0.25; // Always something to eat, but not much.
#ifdef __Hare1950s
          bool veg_patchy = true;
#else   // If it is not the special case of the 1950s
          //
          bool veg_patchy = SupplyVegPatchy(a_polygon);
#endif
          if (veg_patchy)
          {
              // Patchy vegetation - normally full access
              if (veg_height>50)
              {
                  // no food at only at very very tall
                  access -= ((veg_height - 50)* g_VegHeightForageReduction);
                  if (access<0) access = 0;
              }
          }
          else
          {
              if (veg_height>g_FarmIntensivenessH)
              {
                  access -= ((veg_height - g_FarmIntensivenessH)* /* g_FarmIntensiveness * */ g_VegHeightForageReduction);
                  if (access<0) access = 0;
              }
          }
          return access * digest * grazedreduction[SupplyGrazingPressure(a_polygon)];
#endif
}

void Landscape::GISASCII_Output( string outpfile, int UTMX, int UTMY ) {
    FILE* OFILE;
    OFILE = fopen( outpfile.c_str(), "w" );
    if (!OFILE) {
        g_msg->Warn( WARN_FILE, "Landscape::GISASCII_Output() "
            "Unable to open file for writing:",
            outpfile );
        exit( 1 );
    }
    char c = '\n';
    fprintf(OFILE, "ncols         %d\n", m_width);
    fprintf(OFILE, "nrows         %d\n", m_height);
    fprintf(OFILE, "xllcorner     %d\n", UTMX );
    fprintf(OFILE, "yllcorner     %d\n", UTMY );
    fprintf(OFILE, "cellsize      %d\n", 1 );
    fprintf(OFILE, "NODATA_value  %d\n", -9999 );
    // The polyref loop
    if (l_map_ascii_map_entity.value() == 1) {
        for (int y = 0; y < m_height; y++) {
            for (int x = 0; x < m_width; x++) {
                fprintf(OFILE, "%d\t", SupplyPolyRef(x, y));
            }
            fprintf(OFILE, "%c", c );
        }
    }
    // The element type loop
    if (l_map_ascii_map_entity.value() == 2) {
        for (int y = 0; y < m_height; y++) {
            for (int x = 0; x < m_width; x++) {
                fprintf(OFILE, "%d\t", SupplyElementType( x, y ));
            }
            fprintf( OFILE, "%c", c );
        }
    }
    fclose( OFILE );
}

void Landscape::SupplyLEReset(void) {
    le_signal_index = 0;
}

int Landscape::SupplyLENext(void) {
    if ((unsigned int)le_signal_index == m_elems.size()) {
        return -1;
    }
    return m_elems[le_signal_index++]->GetPoly();
}

int Landscape::SupplyLECount(void) {
    return (int)m_elems.size();
}

LE_Signal Landscape::SupplyLESignal(int a_polyref) {
    return m_elems[m_polymapping[a_polyref]]->GetSignal();
}

void Landscape::SetLESignal(int a_polyref, LE_Signal a_signal) {
    m_elems[m_polymapping[a_polyref]]->SetSignal(a_signal);
}


void Landscape::IncTreatCounter(int a_treat) {
    if (a_treat < 0 || a_treat >= last_treatment) {
        char errornum[20];
        sprintf(errornum, "%d", a_treat);
        g_msg->Warn(WARN_BUG, "Landscape::IncTreatCounter(): Index"" out of range!", errornum);
        exit(1);
    }
    m_treatment_counts[a_treat] ++;
}


void Landscape::DumpTreatCounters(const char * a_filename) {
    FILE * l_file = fopen(a_filename, "w");
    if (!l_file) {
        g_msg->Warn(WARN_FILE, "Landscape::DumpTreatCounters(): ""Unable to open file for writing: %s\n", a_filename);
        exit(1);
    }

    for (int i = start; i < last_treatment; i++) {
        fprintf(l_file, "%3d %s %10d\n", i, EventtypeToString(i).c_str(), m_treatment_counts[i]);
    }
    fclose(l_file);
}

void Landscape::DumpMapGraphics(const char * a_filename) {
    unsigned int linesize = m_maxextent * 3;
    unsigned char * frame_buffer = (unsigned char *)malloc(sizeof(unsigned char)* linesize);

    if (frame_buffer == NULL) {
        g_msg->Warn(WARN_FILE, "Landscape::DumpMapGraphics(): Out of memory!", "");
        exit(1);
    }

    FILE * l_file = fopen(a_filename, "w");
    if (!l_file) {
        g_msg->Warn(WARN_FILE, "Landscape::DumpMapGraphics(): ""Unable to open file for writing: %s\n", a_filename);
        exit(1);
    }

    fprintf(l_file, "P6\n%d %d %d\n", m_width, m_height, 255);

    for (int y = 0; y < m_height; y++) {
        int i = 0;
        for (int x = 0; x < m_width; x++) {
            int eletype = (int)SupplyElementType(x, y);
            int localcolor = 16777215 / eletype;

            if (eletype == (int)tole_Field) {
                int category;
                double hei = SupplyVegHeight(x, y);
                if (hei > 50.0) category = 0; else category = (int)(200.0 - (hei * 4.0));
                localcolor = ((category * 65536) + 65535);
            }

            frame_buffer[i++] = (unsigned char)(localcolor & 0xff);
            frame_buffer[i++] = (unsigned char)((localcolor >> 8) & 0xff);
            frame_buffer[i++] = (unsigned char)((localcolor >> 16) & 0xff);
        }
        fwrite(frame_buffer, sizeof(unsigned char), linesize, l_file);
    }

    fclose(l_file);

    free(frame_buffer);
}

void Landscape::FillVegAreaData() {
    for (unsigned int i = 0; i < (tov_Undefined + 1); i++) {
        l_vegtype_areas[i] = 0.0;
    }

    // Sum up statistics on element type.
    for (unsigned int i = 0; i < m_elems.size(); i++) {
        l_vegtype_areas[m_elems[i]->GetVegType()] += m_elems[i]->GetArea();
    }
}

void Landscape::DumpMapInfoByArea(const char * a_filename, bool a_append, bool a_dump_zero_areas, bool a_write_veg_names) {
    FillVegAreaData();
    FILE * outf;
    if (a_append) {
        outf = fopen(a_filename, "a");
        if (!outf) {
            g_msg->Warn(WARN_FILE, "Landscape::DumpMapInfoByArea(): ""Unable to open file for appending", a_filename);
            exit(1);
        }
    }
    else {
        outf = fopen(a_filename, "w");
        if (!outf) {
            g_msg->Warn(WARN_FILE, "Landscape::DumpMapInfoByArea(): ""Unable to open file for writing", a_filename);
            exit(1);
        }
    }

    // Emit element type info.
    for (unsigned int i = 0; i < tov_Undefined + 1; i++) {
        if (i == tov_OFirstYearDanger)
            continue;
        if (!a_dump_zero_areas && l_vegtype_areas[i] < 0.5)
            continue;

        fprintf(outf, "%6ld\t%3d\t%10.0f", g_date->OldDays() + g_date->DayInYear() - 364, i, l_vegtype_areas[i]);
        if (a_write_veg_names)
            fprintf(outf, "\t%s\n", VegtypeToString((TTypesOfVegetation)i).c_str()); else
            fprintf(outf, "\n");
    }

    fclose(outf);
}






LE * Landscape::NewElement(TTypesOfLandscapeElement a_type) {
    LE * elem;
    static char error_num[20];

    switch (a_type) {
    case tole_Hedges:
        elem = new Hedges;
        break;
    case tole_HedgeBank:
        elem = new HedgeBank;
        elem->SetALMaSSEleType(140);
        break;
    case tole_BeetleBank:
        elem = new BeetleBank;
        elem->SetALMaSSEleType(141);
        break;
    case tole_RoadsideVerge:
        elem = new RoadsideVerge;
        break;
    case tole_Railway:
        elem = new Railway;
        break;
    case tole_FieldBoundary:
        elem = new FieldBoundary;
        elem->SetALMaSSEleType(160);
        break;
    case tole_Marsh:
        elem = new Marsh;
        break;
    case tole_Orchard:
        elem = new Orchard;
        break;
    case tole_OrchardBand:
        elem = new OrchardBand;
        break;
    case tole_MownGrass:
        elem = new MownGrass;
        break;
    case tole_Heath:
        elem = new Heath;
        break;
    case tole_Scrub:
        elem = new Scrub;
        break;
    case tole_Field:
        elem = new Field;
        break;
    case tole_PermanentSetaside:
        elem = new PermanentSetaside;
        break;
    case tole_PermPasture:
        elem = new PermPasture;
        break;
    case tole_PermPastureLowYield:
        elem = new PermPastureLowYield;
        break;
    case tole_PermPastureTussocky:
        elem = new PermPastureTussocky;
        break;
    case tole_NaturalGrassDry:
        elem = new NaturalGrassDry;
        break;
    case tole_RiversidePlants:
        elem = new RiversidePlants;
        break;
    case tole_PitDisused:
        elem = new PitDisused;
        break;
    case tole_RiversideTrees:
        elem = new RiversideTrees;
        break;
    case tole_DeciduousForest:
        elem = new DeciduousForest;
        break;
    case tole_MixedForest:
        elem = new MixedForest;
        break;
    case tole_YoungForest:
        elem = new YoungForest;
        break;
    case tole_ConiferousForest:
        elem = new ConiferousForest;
        break;
    case tole_StoneWall:
        elem = new StoneWall;
        break;
    case tole_Garden:
        elem = new Garden;
        break;
    case tole_Track:
        elem = new Track;
        break;
    case tole_SmallRoad:
        elem = new SmallRoad;
        break;
    case tole_LargeRoad:
        elem = new LargeRoad;
        break;
    case tole_Building:
        elem = new Building;
        break;
    case tole_ActivePit:
        elem = new ActivePit;
        break;
    case tole_Freshwater:
        elem = new Freshwater;
        break;
    case tole_Pond:
        elem = new Pond;
        break;
    case tole_River:
        elem = new River;
        break;
    case tole_Saltwater:
        elem = new Saltwater;
        break;
    case tole_Coast:
        elem = new Coast;
        break;
    case tole_BareRock:
        elem = new BareRock;
        break;
    case tole_AmenityGrass:
        elem = new AmenityGrass;
        break;
    case tole_Parkland:
        elem = new Parkland;
        break;
    case tole_UrbanNoVeg:
        elem = new UrbanNoVeg;
        break;
    case tole_UrbanPark:
        elem = new UrbanPark;
        break;
    case tole_BuiltUpWithParkland:
        elem = new BuiltUpWithParkland;
        break;
    case tole_SandDune:
        elem = new SandDune;
        break;
    case tole_Copse:
        elem = new Copse;
        break;
    case tole_RoadsideSlope:
        elem = new RoadsideSlope;
        break;
    case tole_MetalledPath:
        elem = new MetalledPath;
        break;
    case tole_Carpark:
        elem = new Carpark;
        break;
    case tole_Churchyard:
        elem = new Churchyard;
        break;
    case tole_NaturalGrassWet:
        elem = new NaturalGrassWet;
        break;
    case tole_Saltmarsh:
        elem = new Saltmarsh;
        break;
    case tole_Stream:
        elem = new Stream;
        break;
    case tole_HeritageSite:
        elem = new HeritageSite;
        break;
    case tole_UnsprayedFieldMargin:
        elem = new UnsprayedFieldMargin;
        elem->SetALMaSSEleType(31);
        break;
    case tole_Wasteland:
        elem = new Wasteland;
        break;
    case tole_IndividualTree:
        elem = new IndividualTree;
        break;
    case tole_PlantNursery:
        elem = new PlantNursery;
        break;
    case tole_Vildtager:
        elem = new Vildtager;
        break;
    case tole_WindTurbine:
        elem = new WindTurbine;
        break;
    case tole_WoodyEnergyCrop:
        elem = new WoodyEnergyCrop;
        break;
    case tole_WoodlandMargin:
        elem = new WoodlandMargin;
        break;
    case tole_Pylon:
        elem = new Pylon;
        break;
    case tole_FishFarm:
        elem = new FishFarm;
        break;
    case tole_Missing:
        elem = new LE; // These should never be actually used.
        break;
    case tole_Chameleon:
        elem = new ChameleonLE;
        break;
    default:
        sprintf(error_num, "%d", a_type);
        g_msg->Warn(WARN_FILE, "Landscape::NewElement(): Unknown landscape element requested:", error_num);
        exit(1);
    } //switch

    elem->SetElementType(a_type);

    return elem;
}

std::string Landscape::EventtypeToString(int a_event) {
    char error_num[20];

    switch (a_event) {
    case start:
        return "                  start";
    case sleep_all_day:
        return "          sleep_all_day";
    case autumn_plough:
        return "          autumn_plough";
    case autumn_harrow:
        return "          autumn_harrow";
    case autumn_roll:
        return "            autumn_roll";
    case autumn_sow:
        return "             autumn_sow";
    case winter_plough:
        return "          winter_plough";
    case deep_ploughing:
        return "         deep_ploughing";
    case spring_plough:
        return "          spring_plough";
    case spring_harrow:
        return "          spring_harrow";
    case spring_roll:
        return "            spring_roll";
    case spring_sow:
        return "             spring_sow";
    case fp_npks:
        return "                fp_npks";
    case fp_npk:
        return "                 fp_npk";
    case fp_pk:
        return "                  fp_pk";
    case fp_liquidNH3:
        return "           fp_liquidNH3";
    case fp_slurry:
        return "              fp_slurry";
    case fp_manganesesulphate:
        return "   fp_manganesesulphate";
    case fp_manure:
        return "              fp_manure";
    case fp_greenmanure:
        return "         fp_greenmanure";
    case fp_sludge:
        return "              fp_sludge";
    case fa_npk:
        return "                 fa_npk";
    case fa_pk:
        return "                  fa_pk";
    case fa_slurry:
        return "              fa_slurry";
    case fa_ammoniumsulphate:
        return "    fa_ammoniumsulphate";
    case fa_manure:
        return "              fa_manure";
    case fa_greenmanure:
        return "         fa_greenmanure";
    case fa_sludge:
        return "              fa_sludge";
    case herbicide_treat:
        return "        herbicide_treat";
    case growth_regulator:
        return "       growth_regulator";
    case fungicide_treat:
        return "        fungicide_treat";
    case insecticide_treat:
        return "      insecticide_treat";
    case product_treat:
        return "pesticide_product_treat";
    case syninsecticide_treat:
        return "   syninsecticide_treat";
    case molluscicide:
        return "           molluscicide";
    case row_cultivation:
        return "        row_cultivation";
    case strigling:
        return "              strigling";
    case flammebehandling:
        return "       flammebehandling";
    case hilling_up:
        return "             hilling_up";
    case water:
        return "                  water";
    case swathing:
        return "               swathing";
    case harvest:
        return "                harvest";
    case cattle_out:
        return "             cattle_out";
    case pigs_out:
        return "               pigs_out";
    case cut_to_hay:
        return "             cut_to_hay";
    case cut_to_silage:
        return "          cut_to_silage";
    case straw_chopping:
        return "         straw_chopping";
    case hay_turning:
        return "            hay_turning";
    case hay_bailing:
        return "            hay_bailing";
    case stubble_harrowing:
        return "      stubble_harrowing";
    case autumn_or_spring_plough:
        return "autumn_or_spring_plough";
    case burn_straw_stubble:
        return "     burn_straw_stubble";
    case mow:
        return "                    mow";
    case cut_weeds:
        return "              cut_weeds";
    case strigling_sow:
        return "          strigling_sow";
    case trial_insecticidetreat:
        return "PesticideTrialTreatment";
    case trial_toxiccontrol:
        return "    PesticideTrialToxic";
    case trial_control:
        return "  PesticideTrialControl";
    case glyphosate:
        return " Glyphosate on setaside";
    default:
        sprintf(error_num, "%d", a_event);
        g_msg->Warn(WARN_FILE, "Landscape::EventtypeToString(): Unknown event type:", error_num);
        exit(1);
    }
}



std::string Landscape::PolytypeToString(TTypesOfLandscapeElement a_le_type) {
    char error_num[20];

    switch (a_le_type) {
    case tole_Hedges:
        return "                Hedge";
    case tole_RoadsideVerge:
        return "       Roadside Verge";
    case tole_Railway:
        return "              Railway";
    case tole_FieldBoundary:
        return "       Field Boundary";
    case tole_Marsh:
        return "                Marsh";
    case tole_Scrub:
        return "                Scrub";
    case tole_Field:
        return "                Field";
    case tole_PermPastureTussocky:
        return " PermPastureTussocky";
    case tole_PermanentSetaside:
        return "   Permanent Setaside";
    case tole_PermPasture:
        return "    Permanent Pasture";
    case tole_PermPastureLowYield:
        return " PermPastureLowYield";
    case tole_NaturalGrassDry:
        return "        Natural Grass";
    case tole_NaturalGrassWet:
        return "    Natural Grass Wet";
    case tole_RiversidePlants:
        return "     Riverside Plants";
    case tole_PitDisused:
        return "          Pit Disused";
    case tole_RiversideTrees:
        return "      Riverside Trees";
    case tole_DeciduousForest:
        return "     Deciduous Forest";
    case tole_MixedForest:
        return "         Mixed Forest";
    case tole_ConiferousForest:
        return "    Coniferous Forest";
    case tole_YoungForest:
        return "         Young Forest";
    case tole_StoneWall:
        return "           Stone Wall";
    case tole_Garden:
        return "               Garden";
    case tole_Track:
        return "                Track";
    case tole_SmallRoad:
        return "           Small Road";
    case tole_LargeRoad:
        return "           Large Road";
    case tole_Building:
        return "             Building";
    case tole_ActivePit:
        return "           Active Pit";
    case tole_Pond:
        return "                 Pond";
    case tole_Freshwater:
        return "          Fresh Water";
    case tole_River:
        return "                River";
    case tole_Saltwater:
        return "            Saltwater";
    case tole_Coast:
        return "                Coast";
    case tole_BareRock:
        return "            Bare Rock";
    case tole_HedgeBank:
        return "            Hedgebank";
    case tole_Heath:
        return "                Heath";
    case tole_Orchard:
        return "              Orchard";
    case tole_OrchardBand:
        return "         Orchard Band";
    case tole_MownGrass:
        return "           Mown Grass";
    case tole_UnsprayedFieldMargin:
        return " UnsprayedFieldMargin";
    case tole_AmenityGrass:
        return "         AmenityGrass";
    case tole_Parkland:
        return "             Parkland";
    case tole_UrbanNoVeg:
        return "           UrbanNoVeg";
    case tole_UrbanPark:
        return "            UrbanPark";
    case tole_BuiltUpWithParkland:
        return "  BuiltUpWithParkland";
    case tole_SandDune:
        return "             SandDune";
    case tole_Copse:
        return "                Copse";
    case tole_RoadsideSlope:
        return "        RoadsideSlope";
    case tole_MetalledPath:
        return "         MetalledPath";
    case tole_Carpark:
        return "             Carpark";
    case tole_Churchyard:
        return "          Churchyard";
    case tole_Saltmarsh:
        return "           Saltmarsh";
    case tole_Stream:
        return "              Stream";
    case tole_HeritageSite:
        return "        HeritageSite";
    case tole_BeetleBank:  //141
        return "         Beetle Bank";
    case tole_UnknownGrass:
        return "       Unknown Grass";
    case tole_Wasteland:
        return " Waste/Building Land";
    case tole_IndividualTree:
        return "      IndividualTree";
    case tole_PlantNursery:
        return "        PlantNursery";
    case tole_Vildtager:
        return "           Vildtager";
    case tole_WindTurbine:
        return "         WindTurbine";
    case tole_WoodyEnergyCrop:
        return "     WoodyEnergyCrop";
    case tole_WoodlandMargin:
        return "      WoodlandMargin";
    case tole_Pylon:
        return "               Pylon";
    case tole_FishFarm:
        return "            FishFarm";
    case tole_Foobar:
    default:
        sprintf(error_num, "%d", a_le_type);
        g_msg->Warn(WARN_FILE, "Landscape::PolytypeToString(): Unknown event type:", error_num);
        exit(1);
    }
}



std::string Landscape::VegtypeToString(TTypesOfVegetation a_veg) {
    char error_num[20];

    switch (a_veg) {
    case tov_Carrots:
        return "Carrots              ";
    case tov_BroadBeans:
        return "BroadBeans           ";
    case tov_FodderGrass:
        return "FodderGrass          ";
    case tov_CloverGrassGrazed1:
        return "CloverGrassGrazed1   ";
    case tov_CloverGrassGrazed2:
        return "CloverGrassGrazed2   ";
    case tov_FieldPeas:
        return "FieldPeas            ";
    case tov_FieldPeasSilage:
        return "FieldPeasSilage      ";
    case tov_FodderBeet:
        return "FodderBeet           ";
    case tov_SugarBeet:
        return "SugarBeet            ";
    case tov_OFodderBeet:
        return "OFodderBeet          ";
    case tov_Lawn:
        return "Lawn                 ";
    case tov_Maize:
        return "Maize                ";
    case tov_MaizeSilage:
        return "MaizeSilage          ";
    case tov_OMaizeSilage:
        return "OMaizeSilage         ";
    case tov_NaturalGrass:
        return "NaturalGrass         ";
    case tov_NoGrowth:
        return "NoGrowth             ";
    case tov_None:
        return "None                 ";
    case tov_OrchardCrop:
        return "OrchardCrop          ";
    case tov_Oats:
        return "Oats                 ";
    case tov_OBarleyPeaCloverGrass:
        return "OBarleyPeaCloverGrass";
    case tov_OCarrots:
        return "OCarrots             ";
    case tov_OCloverGrassGrazed1:
        return "OCloverGrassGrazed1  ";
    case tov_OCloverGrassGrazed2:
        return "OCloverGrassGrazed2  ";
    case tov_OCloverGrassSilage1:
        return "OCloverGrassSilage1  ";
    case tov_OFieldPeas:
        return "OFieldPeas           ";
    case tov_OFieldPeasSilage:
        return "OFieldPeasSilage     ";
    case tov_OFirstYearDanger:
        return "OFirstYearDanger     ";
    case tov_OGrazingPigs:
        return "OGrazingPigs         ";
    case tov_OOats:
        return "OOats                ";
    case tov_OPermanentGrassGrazed:
        return "OPermanentGrassGrazed";
    case tov_OPotatoes:
        return "OPotatoes            ";
    case tov_OSBarleySilage:
        return "OSBarleySilage       ";
    case tov_OSeedGrass1:
        return "OSeedGrass1          ";
    case tov_OSeedGrass2:
        return "OSeedGrass2          ";
    case tov_OSetaside:
        return "OSetaside             ";
    case tov_OSpringBarley:
        return "OSpringBarley        ";
    case tov_OSpringBarleyExt:
        return "OSpringBarleyExt     ";
    case tov_OSpringBarleyClover:
        return "OSpringBarleyClover  ";
    case tov_OSpringBarleyGrass:
        return "OSpringBarleyGrass   ";
    case tov_OSpringBarleyPigs:
        return "OSpringBarleyPigs    ";
    case tov_OTriticale:
        return "OTriticale           ";
    case tov_OWinterBarley:
        return "OWinterBarley        ";
    case tov_OWinterBarleyExt:
        return "OWinterBarleyExt     ";
    case tov_OWinterRape:
        return "OWinterRape          ";
    case tov_OWinterRye:
        return "OWinterRye           ";
    case tov_OWinterWheatUndersown:
        return "OWinterWheatUndersown";
    case tov_OWinterWheat:
        return "OWinterWheat";
    case tov_OWinterWheatUndersownExt:
        return "OWinterWheatUsowExt  ";
    case tov_PermanentGrassGrazed:
        return "PermanentGrassGrazed ";
    case tov_PermanentGrassLowYield:
        return "PermanentGrassLowYield";
    case tov_PermanentGrassTussocky:
        return "PermanentGrassTussocky";
    case tov_PermanentSetaside:
        return "PermanentSetaside    ";
    case tov_Potatoes:
        return "PotatoesEat          ";
    case tov_PotatoesIndustry:
        return "PotatoesIndustry     ";
    case tov_SeedGrass1:
        return "SeedGrass1           ";
    case tov_SeedGrass2:
        return "SeedGrass2           ";
    case tov_Setaside:
        return "Setaside             ";
    case tov_SpringBarley:
        return "SpringBarley         ";
    case tov_SpringBarleySpr:
        return "SpringBarleySpr      ";
    case tov_SpringBarleyPTreatment:
        return "SpringBarleyPTreat   ";
    case tov_SpringBarleySKManagement:
        return "SpringBarleySKMan    ";
    case tov_SpringBarleyCloverGrass:
        return "SprBarleyCloverGrass ";
    case tov_SpringBarleyGrass:
        return "SpringBarleyGrass    ";
    case tov_SpringBarleySeed:
        return "SpringBarleySeed     ";
    case tov_SpringBarleySilage:
        return "SpringBarleySilage   ";
    case tov_SpringRape:
        return "SpringRape           ";
    case tov_SpringWheat:
        return "SpringWheat          ";
    case tov_AgroChemIndustryCereal:
        return "AgroChemIndustry Cereal      ";
    case tov_Triticale:
        return "Triticale            ";
    case tov_WinterBarley:
        return "WinterBarley         ";
    case tov_WinterRape:
        return "WinterRape           ";
    case tov_WinterRye:
        return "WinterRye            ";
    case tov_WinterWheat:
        return "WinterWheat          ";
    case tov_WinterWheatShort:
        return "WinterWheatShort     ";
    case tov_WWheatPControl:
        return "P Trial Control      ";
    case tov_WWheatPToxicControl:
        return "P Trial Toxic Control";
    case tov_WWheatPTreatment:
        return "P Trial Treatment    ";
    case tov_Undefined:
        return "Undefined            ";
    case tov_WinterWheatStrigling:
        return "WWStrigling          ";
    case tov_WinterWheatStriglingSingle:
        return "WWStriglingSingle    ";
    case tov_WinterWheatStriglingCulm:
        return "WWStriglingCulm      ";
    case tov_SpringBarleyCloverGrassStrigling:
        return "SBPCGStrigling       ";
    case tov_SpringBarleyStrigling:
        return "SBarleyStrigling     ";
    case tov_SpringBarleyStriglingSingle:
        return "SBarleyStriglingSgl  ";
    case tov_SpringBarleyStriglingCulm:
        return "SBarleyStriglingCulm ";
    case tov_MaizeStrigling:
        return "MaizseStrigling      ";
    case tov_WinterRapeStrigling:
        return "WRapeStrigling       ";
    case tov_WinterRyeStrigling:
        return "WRyeStrigling        ";
    case tov_WinterBarleyStrigling:
        return "WBStrigling          ";
    case tov_FieldPeasStrigling:
        return "FieldPeasStrigling   ";
    case tov_SpringBarleyPeaCloverGrassStrigling:
        return "SBPeaCloverGrassStr  ";
    case tov_YoungForest:
        return "Young Forest         ";
    case tov_Wasteland:
        return "Wasteland            ";
    case tov_Heath:
        return "Heath/Grass          ";
    case  tov_PlantNursery:
        return "Plant Nursery        ";

    default:
        sprintf(error_num, "%d", a_veg);
        g_msg->Warn(WARN_FILE, "Landscape::VegtypeToString(): Unknown event type:", error_num);
        exit(1);
    }
}



bool Landscape::BorderNeed(TTypesOfLandscapeElement a_letype) {
    static char error_num[20];
    bool AddBorder = false;
    switch (a_letype) {
        // No border is needed toward these neighbouring element types.
    case tole_Hedges:
    case tole_HedgeBank:
    case tole_BeetleBank:
    case tole_RoadsideVerge:
    case tole_Marsh:
    case tole_RiversidePlants:
    case tole_UnsprayedFieldMargin:
    case tole_OrchardBand:
    case tole_MownGrass:
        break;

    case tole_IndividualTree:
    case tole_PlantNursery:
    case tole_Vildtager:
    case tole_WindTurbine:
    case tole_WoodyEnergyCrop:
    case tole_WoodlandMargin:
    case tole_Pylon:
    case tole_NaturalGrassDry:
    case tole_Railway:
    case tole_FieldBoundary:
    case tole_Scrub:
    case tole_Field:
    case tole_PermanentSetaside:
    case tole_PermPasture:
    case tole_PermPastureTussocky:
    case tole_PermPastureLowYield:
    case tole_PitDisused:
    case tole_RiversideTrees:
    case tole_DeciduousForest:
    case tole_MixedForest:
    case tole_YoungForest:
    case tole_ConiferousForest:
    case tole_StoneWall:
    case tole_Garden:
    case tole_Track:
    case tole_SmallRoad:
    case tole_LargeRoad:
    case tole_Building:
    case tole_ActivePit:
    case tole_Pond:
    case tole_FishFarm:
    case tole_Freshwater:
    case tole_River:
    case tole_Saltwater:
    case tole_Coast:
    case tole_BareRock:
    case tole_Heath:
    case tole_Orchard:
    case tole_AmenityGrass:
    case tole_Parkland:
    case tole_UrbanNoVeg:
    case tole_UrbanPark:
    case tole_BuiltUpWithParkland:
    case tole_SandDune:
    case tole_Copse:
    case tole_NaturalGrassWet:
    case tole_RoadsideSlope:
    case tole_MetalledPath:
    case tole_Carpark:
    case tole_Churchyard:
    case tole_Saltmarsh:
    case tole_Stream:
    case tole_HeritageSite:
        AddBorder = true;
        break;

    default:
        sprintf(error_num, "%d", a_letype);
        g_msg->Warn(WARN_BUG, "Landscape::BorderNeed(): Unknown element type:", error_num);
        exit(1);
    }
    return AddBorder;
}



/* Caution: The trace file file generated has an offset * in the reported date off by one! This looks like
* something is severely broken, but in reality I'm just * too lazy to fix this (design problem). */
/*
void Landscape::TestCropManagement( void ) {
FILE * tracefile;

tracefile = fopen(g_farm_test_crop_filename.value(), "w" )
if (!tracefile){
g_msg->Warn( WARN_FILE, "Landscape::TestCropManagement(): ""Unable to open file for writing: ",
g_farm_test_crop_filename.value() );
exit( 1 );
}

// Loop through the polygons and find one with 'our'
// vegetation type.
int l_poly;
for ( unsigned int i = 0; i < m_elems.size(); i++ ) {
if ( m_elems[ i ]->GetVegType() == g_farm_test_crop_type.value() ) {
l_poly = m_elems[ i ]->GetPoly();
break;
}
}

//int poly = SupplyPolyRef( g_farm_test_crop_sample_x.value(),
//              g_farm_test_crop_sample_y.value() );
//int poly = 943;
int inf = 0;
int d = 0;
while ( 1 ) {  //Mysterious - what does this do?
Tick();
if ( ++inf > 364 ) {
printf( "%6d\n", ++d );
inf = 0;
}
}

fprintf( tracefile, "Polygon: %d\n", l_poly );
fprintf( tracefile, "EleType: %d\n", SupplyElementType( l_poly ) ),
fprintf( tracefile, "VegType: %d\n", SupplyVegType( l_poly ) );

for ( int x = 0; x < g_farm_test_crop_daystorun.value(); x++ ) {
Tick();
int index = 0, last;
int day = SupplyDayInMonth();
int month = SupplyMonth();
int year = SupplyYear();
long date = SupplyGlobalDate();
int veg = SupplyVegType( l_poly );

do {
last = SupplyLastTreatment( l_poly, & index );

// Comment out if you want to know *everything*,
// including nothing, that happens evey day on the test field.
if ( last == sleep_all_day && !g_farm_test_crop_reportdaily.value() )
continue;

fprintf( tracefile, "%2d %2d %4d %5ld %s %2d :: ""%7.2f %7.2f %7.2f %7.2f\n", day, month, year, date, EventtypeToString( last ),
veg, SupplyVegBiomass( l_poly ), SupplyVegCover( l_poly ), SupplyVegHeight( l_poly ), SupplyInsects( l_poly ) );
} while ( last != sleep_all_day );
}

fclose( tracefile );
}
*/

void Landscape::DumpVegAreaData(int a_day) {

    if (cfg_dumpvegjan.value()) {
        if ((a_day % 365) == 0) { // Jan 1st
            DumpMapInfoByArea(cfg_dumpvegjanfile.value(), true, true, true);
            return;
        }
    }
    if (cfg_dumpvegjune.value()) {
        if ((a_day % 365) == 152) { // 1st June
            DumpMapInfoByArea(cfg_dumpvegjunefile.value(), true, true, true);
        }
    }

}
//------------------------------------------------------------------------------