Creating small area income and small area coastal income maps for Ireland, Northern Ireland, Scotland, Wales, England, France, and Belgium

CoastalIncomeAnalysis.cs

using CsvHelper;
using GeoAPI.CoordinateSystems.Transformations;
using NetTopologySuite.Features;
using NetTopologySuite.Geometries;
using NetTopologySuite.Geometries.Prepared;
using NetTopologySuite.IO;
using NetTopologySuite.Operation.Buffer;
using Newtonsoft.Json;
using NUnit.Framework;
using ProjNet.CoordinateSystems;
using ProjNet.CoordinateSystems.Transformations;
using System;
using System.Collections.Generic;
using System.Globalization;
using System.IO;
using System.Linq;
using System.Text;

namespace CoastalEconomies
{
    class Program
    {
        static void Main(string[] args)
        {
            // Producing Shapefiles with NetTopologySuite
            // https://gist.github.com/billpratt/6bf4ab8e4435dc39db92

            // You can learn more about projections at epsg.io
            //string EPSG3035WKT = File.ReadAllText("ProjectionWKT/EPSG3035.txt"); // this is the best North Europe projection with units in metres, but it isn't supported by ProjNet
            string EPSG23030WKT = File.ReadAllText("ProjectionWKT/EPSG23030.txt"); // this is a pretty good North Europe projection with units in metres that is supported by ProjNet
                                                                                   //string EPSG3857WKT = File.ReadAllText("ProjectionWKT/EPSG3857.txt"); // this is WGS84 (not sure what more)
                                                                                   //string EPSG4326WKT = File.ReadAllText("ProjectionWKT/EPSG4326.txt"); // this is WGS84 (not sure what more)

            CoordinateSystemFactory csf = new CoordinateSystemFactory();
            CoordinateTransformationFactory trf = new CoordinateTransformationFactory();
            ICoordinateTransformation TransformToMetres = trf.CreateFromCoordinateSystems(GeographicCoordinateSystem.WGS84, csf.CreateFromWkt(EPSG23030WKT));
            ICoordinateTransformation TransformToDegrees = trf.CreateFromCoordinateSystems(csf.CreateFromWkt(EPSG23030WKT), GeographicCoordinateSystem.WGS84);

            // re-usable geojsonreader
            GeoJsonReader gjsonreader = new NetTopologySuite.IO.GeoJsonReader();

            Console.WriteLine("Load France commune boundaries.");
            FeatureCollection FranceCommuneBoundaries = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../FranceCommunes_4326.geojson"));

            Console.WriteLine("Load Scotland datazone boundaries.");
            FeatureCollection ScotlandDataZoneBoundaries = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../ScotlandDatazones_4326.geojson"));

            Console.WriteLine("Load England and Wales MSOA boundaries.");
            FeatureCollection EandWMSOABoundaries = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../EandWMSOAs_4326.geojson"));

            Console.WriteLine("Load Belgium statistical sector boundaries.");
            FeatureCollection BelgiumSSBoundaries = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../sh_statbel_statistical_sectors_4326.geojson"));

            Console.WriteLine("Load Northern Ireland super output area boundaries.");
            FeatureCollection NISOAs = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../NISOAs_4326.geojson"));

            Console.WriteLine("Load Republic of Ireland electoral division boundaries."); // http://census2016.geohive.ie/search?tags=Census%202016%20Boundaries
            FeatureCollection ROIEDs = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../ROI_ElectoralDivision_4326.geojson"));

            // Belgian data is at https://twitter.com/thomasforth/status/1229478951682269186

            //Console.WriteLine("Load Ireland coastline.");
            FeatureCollection IrelandCoastFeatures = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../IrelandCoast_4326.geojson"));

            // load coastline of GBandICoast
            //FeatureCollection GBandICoast = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../GBAndIreland_4326.geojson"));

            Console.WriteLine("Load North West Europe (simplified) coastline.");
            FeatureCollection NWEuropeCoastFeatures = gjsonreader.Read<FeatureCollection>(File.ReadAllText(@"../../../../../FullCoastAtFullQuality.geojson"));

            // Create a rectangle for clipping the coastline (the projection we use is only valid in North West Europe).
            Console.WriteLine("Clipping the loaded coastline to area of interest (and safe to use with EPSG23030 projection).");
            List<Coordinate> coords = new List<Coordinate>();

            // -11.53,40.69,10.36,61.0
            coords.Add(new Coordinate() { X = -12, Y = 61 });
            coords.Add(new Coordinate() { X = 11, Y = 61 });
            coords.Add(new Coordinate() { X = 11, Y = 40 });
            coords.Add(new Coordinate() { X = -12, Y = 40 });
            coords.Add(new Coordinate() { X = -12, Y = 61 });

            Polygon pg = new Polygon(new LinearRing(coords.ToArray()));

            FeatureCollection BufferedCoastCollection = new FeatureCollection();
            Console.WriteLine("Buffering and dissolving the coastline.");
            FeatureCollection InvalidCoastlineFeatures = new FeatureCollection();
            foreach (Feature coast in NWEuropeCoastFeatures)
            {
                try
                {
                    coast.Geometry = pg.Intersection(coast.Geometry.Boundary);
                    Geometry coastinmetres = Transform(coast.Geometry, (MathTransform)TransformToMetres.MathTransform);
                    Geometry bufferedcoast = new BufferOp(coastinmetres).GetResultGeometry(500);
                    bufferedcoast = new BufferOp(bufferedcoast).GetResultGeometry(9500);
                    Geometry bufferedcoastindegrees = Transform(bufferedcoast, (MathTransform)TransformToDegrees.MathTransform);
                    if (bufferedcoastindegrees.IsValid)
                    {
                        BufferedCoastCollection.Add(new Feature(bufferedcoastindegrees, null));
                    }
                    else
                    {
                        InvalidCoastlineFeatures.Add(coast);
                    }
                }
                catch
                {
                    InvalidCoastlineFeatures.Add(coast);
                }
            }

            GeoJsonWriter writer = new GeoJsonWriter();
            writer.SerializerSettings.Formatting = Formatting.Indented;

            // The GeoJSON standard says you should WGS84 and nothing else as a projection -- https://tools.ietf.org/html/rfc7946#section-4

            File.WriteAllText($"InvalidCoastlineFeatures.geojson", writer.Write(InvalidCoastlineFeatures));

            //File.WriteAllText($"BufferedCoastDebug.geojson", writer.Write(BufferedCoastCollection));

            Console.WriteLine("Writing buffered and coastal area for debug.");
            Geometry DissolvedBufferedCoast = NetTopologySuite.Operation.Union.CascadedPolygonUnion.Union(new List<Geometry>(BufferedCoastCollection.Select(x => x.Geometry)));
            File.WriteAllText($"BufferedCoastDebug_unioned.geojson", writer.Write(DissolvedBufferedCoast));

            PreparedPolygon IndexedCoast = new PreparedPolygon((IPolygonal)DissolvedBufferedCoast);
            NWEuropeCoastFeatures.Clear();
            BufferedCoastCollection.Clear();

            Console.WriteLine("Loading Republic of Ireland Electoral District income data.");
            List<ROIIncome> ROIIncomes = new List<ROIIncome>();
            using (StreamReader reader = new StreamReader(@"Data/IIA01.csv"))
            {
                using (CsvReader csv = new CsvReader(reader, CultureInfo.InvariantCulture))
                {
                    ROIIncomes = csv.GetRecords<ROIIncome>().ToList();
                }
            }


            Console.WriteLine("Loading Northern Ireland super output area income deprivation index.");
            Dictionary<string, NIIncome> NIIncomes = new Dictionary<string, NIIncome>();
            using (StreamReader reader = new StreamReader(@"Data/NISOA_income_Deprivation.csv"))
            {
                using (CsvReader csv = new CsvReader(reader, CultureInfo.InvariantCulture))
                {
                    NIIncomes = csv.GetRecords<NIIncome>().ToDictionary(x => x.SOA2001, x => x);
                }
            }

            
            Console.WriteLine("Loading Scotland datazone income data.");
            Dictionary<string, ScotlandIncome> ScotlandIncomes = new Dictionary<string, ScotlandIncome>();
            using (StreamReader reader = new StreamReader(@"Data/ScotlandDZIncomes.csv"))
            {
                using (CsvReader csv = new CsvReader(reader, CultureInfo.InvariantCulture))
                {
                    ScotlandIncomes = csv.GetRecords<ScotlandIncome>().ToDictionary(x => x.S2011DataZonecode, x => x);
                }
            }

            Console.WriteLine("Loading England and Wales MSOA income data.");
            Dictionary<string, EandWIncome> EandWIncomes = new Dictionary<string, EandWIncome>();
            using (StreamReader reader = new StreamReader(@"Data/EandWIncome.csv"))
            {
                using (CsvReader csv = new CsvReader(reader, CultureInfo.InvariantCulture))
                {
                    EandWIncomes = csv.GetRecords<EandWIncome>().ToDictionary(x => x.MSOAcode, x => x);
                }
            }

            Console.WriteLine("Loading France commune income data.");
            Dictionary<string, FranceIncome> FranceIncomes = new Dictionary<string, FranceIncome>();
            using (StreamReader reader = new StreamReader(@"Data/cc_filosofi_2017_COM.CSV"))
            {
                using (CsvReader csv = new CsvReader(reader, CultureInfo.InvariantCulture))
                {
                    csv.Configuration.Delimiter = ";";
                    FranceIncomes = csv.GetRecords<FranceIncome>().ToDictionary(x => x.CODGEO, x => x);
                }
            }

            Console.WriteLine("Loading Belgium statistical sector income data.");
            Dictionary<string, BelgiumIncome> BelgIncomes = new Dictionary<string, BelgiumIncome>();
            using (StreamReader reader = new StreamReader(@"Data/BelgIncome.csv"))
            {
                using (CsvReader csv = new CsvReader(reader, CultureInfo.InvariantCulture))
                {
                    BelgIncomes = csv.GetRecords<BelgiumIncome>().ToDictionary(x => x.Code, x => x);
                }
            }

            FeatureCollection AllGeographies = new FeatureCollection();
            FeatureCollection InvalidPolygons = new FeatureCollection();

            Console.WriteLine("Calculating Northern Ireland SOAs that overlap the 'coastal' area.");
            foreach (Feature nisoa in NISOAs)
            {
                string soacode = (string)nisoa.Attributes["SOA_CODE"];
                nisoa.Attributes.GetNames().ToList().ForEach(x => nisoa.Attributes.DeleteAttribute(x));
                
                nisoa.Attributes.Add("id", soacode);
                nisoa.Attributes.Add("country", "Northern Ireland");
                nisoa.Attributes.Add("income_mean", 0);
                nisoa.Attributes.Add("income_mdn", 0);
                nisoa.Attributes.Add("income_rtn", 0.0);
                nisoa.Attributes.Add("coastal", false);

                Geometry TestGeometry = nisoa.Geometry;
                if (TestGeometry.IsValid == false)
                {
                    TestGeometry = new BufferOp(nisoa.Geometry).GetResultGeometry(0);
                    Assert.IsTrue(TestGeometry.IsValid);
                }
                if (NIIncomes.ContainsKey(soacode))
                {
                    nisoa.Attributes["income_rtn"] = (double)NIIncomes[soacode].IncomeDomainRank_1mostdeprived;
                }
                if (IndexedCoast.Intersects(TestGeometry))
                {
                    nisoa.Attributes["coastal"] = true;
                    //CoastalGeographies.Add(nisoa);
                }
                Assert.IsTrue(nisoa.Attributes.GetType("income_rtn") == typeof(double));
                AllGeographies.Add(nisoa);
            }

            Console.WriteLine("Calculating Republic of Ireland electoral districts that overlap the 'coastal' area.");
            double ROIAverageIncome = ROIIncomes.Average(x => x.MedianIncome);
            foreach (Feature roied in ROIEDs)
            {
                string name = (string)roied.Attributes["CSOED_34_1"];
                string id = (string)roied.Attributes["CSOED_3409"];

                roied.Attributes.GetNames().ToList().ForEach(x => roied.Attributes.DeleteAttribute(x));
                roied.Attributes.Add("id", id + name);
                roied.Attributes.Add("country", "Ireland");
                roied.Attributes.Add("income_mean", 0);
                roied.Attributes.Add("income_mdn", 0);
                roied.Attributes.Add("income_rtn", 0.0);
                roied.Attributes.Add("coastal", false);

                id = id.Substring(2);

                ROIIncome income;
                if (id.Contains("/"))
                {
                    income = ROIIncomes.Where(x => x.ED.Contains(name)).FirstOrDefault();
                }
                else
                {
                    income = ROIIncomes.Where(x => x.ED.StartsWith(id) && x.ED.Contains(name)).FirstOrDefault();
                }

                if (income != null)
                {
                    roied.Attributes["income_rtn"] = (double)(100 * income.MedianIncome / ROIAverageIncome);
                    roied.Attributes["income_mdn"] = income.MedianIncome;
                }

                if (IndexedCoast.Intersects(roied.Geometry))
                {
                    roied.Attributes["coastal"] = true;
                    //CoastalGeographies.Add(roied);
                }
                Assert.IsTrue(roied.Attributes.GetType("income_rtn") == typeof(double));
                AllGeographies.Add(roied);
            }

            
            Console.WriteLine("Calculating England and Wales MSOAs that overlap the 'coastal' area.");
            double EandWAverageIncome = EandWIncomes.Average(x => x.Value.NAIBHC);
            foreach (Feature msoa in EandWMSOABoundaries)
            {
                Geometry TestGeometry = msoa.Geometry;
                if (TestGeometry.IsValid == false)
                {
                    TestGeometry = new BufferOp(msoa.Geometry).GetResultGeometry(0);
                    Assert.IsTrue(TestGeometry.IsValid);
                }

                string msoacode = (string)msoa.Attributes["msoa11cd"];
                double income = EandWIncomes[msoacode].NAIBHC;                
                msoa.Attributes.GetNames().ToList().ForEach(x => msoa.Attributes.DeleteAttribute(x));

                msoa.Attributes.Add("id", msoacode);
                if (msoacode.StartsWith("W"))
                {
                    msoa.Attributes.Add("country", "Wales");
                }
                else
                {
                    msoa.Attributes.Add("country", "England");
                }
                msoa.Attributes.Add("income_rtn", (double)(100 * income / EandWAverageIncome));
                msoa.Attributes.Add("income_mean", 0);
                msoa.Attributes["income_mean"] = (int)income;
                msoa.Attributes.Add("income_mdn", 0);
                msoa.Attributes.Add("coastal", false);

                

                if (IndexedCoast.Intersects(TestGeometry))
                {
                    msoa.Attributes["coastal"] = true;
                    //CoastalGeographies.Add(msoa);
                }
                Assert.IsTrue(msoa.Attributes.GetType("income_rtn") == typeof(double));
                AllGeographies.Add(msoa);
            }

            Console.WriteLine("Calculating Belgium statistical sectors that overlap the 'coastal' area.");
            BelgIncomes = new Dictionary<string, BelgiumIncome>(BelgIncomes.Where(x => x.Value.Revenumedianpardec.HasValue));
            double BelgiumAverageIncome = BelgIncomes.Average(x => x.Value.Revenumedianpardec.Value);
            foreach (Feature ss in BelgiumSSBoundaries)
            {
                Geometry TestGeometry = ss.Geometry;
                if (TestGeometry.IsValid == false)
                {
                    TestGeometry = new BufferOp(ss.Geometry).GetResultGeometry(0);
                    Assert.IsTrue(TestGeometry.IsValid);
                }
                string sectorcode = (string)ss.Attributes["CD_SECTOR"];
                ss.Attributes.GetNames().ToList().ForEach(x => ss.Attributes.DeleteAttribute(x));

                ss.Attributes.Add("id", sectorcode);
                ss.Attributes.Add("country", "Belgium");
                ss.Attributes.Add("income_rtn", 0.0);
                ss.Attributes.Add("income_mean", 0);
                ss.Attributes.Add("income_mdn", 0);
                ss.Attributes.Add("coastal", false);

                if (BelgIncomes.ContainsKey(sectorcode))
                {
                    double income = BelgIncomes[sectorcode].Revenumedianpardec.Value;
                    ss.Attributes["income_mdn"] = BelgIncomes[sectorcode].Revenumedianpardec.Value;
                    if (BelgIncomes[sectorcode].Revenumoyenpardec.HasValue)
                    {
                        ss.Attributes["income_mean"] = BelgIncomes[sectorcode].Revenumoyenpardec.Value;
                    }
                    ss.Attributes["income_rtn"] = (double)(100 * income / BelgiumAverageIncome);
                    Assert.IsTrue(ss.Attributes.GetType("income_rtn") == typeof(double));

                    if (IndexedCoast.Intersects(TestGeometry))
                    {
                        ss.Attributes["coastal"] = true;
                        //CoastalGeographies.Add(ss);
                    }
                    AllGeographies.Add(ss);
                }
            }


            Console.WriteLine("Calculating Scottish datazones that overlap the 'coastal' area.");
            double ScotlandAverageIncome = ScotlandIncomes.Average(x => x.Value.MedianGHIWeek);
            foreach (Feature SDZ in ScotlandDataZoneBoundaries)
            {
                string datazonecode = (string)SDZ.Attributes["DataZone"];
                SDZ.Attributes.GetNames().ToList().ForEach(x => SDZ.Attributes.DeleteAttribute(x));
                SDZ.Attributes.Add("id", datazonecode);
                SDZ.Attributes.Add("country", "Scotland");
                SDZ.Attributes.Add("income_rtn", 0.0);
                SDZ.Attributes.Add("income_mean", 0);
                SDZ.Attributes.Add("income_mdn", 0);
                SDZ.Attributes.Add("coastal", false);

                Geometry TestGeometry = SDZ.Geometry;
                if (TestGeometry.IsValid == false)
                {
                    TestGeometry = new BufferOp(SDZ.Geometry).GetResultGeometry(0);
                    if (TestGeometry.IsValid == false)
                    {
                        InvalidPolygons.Add(SDZ);
                    }
                }

                if (ScotlandIncomes.ContainsKey(datazonecode))
                {
                    double income = ScotlandIncomes[datazonecode].MedianGHIWeek;
                    SDZ.Attributes["income_rtn"] = (double)(100 * income / ScotlandAverageIncome);
                    SDZ.Attributes["income_mdn"] = 52* ScotlandIncomes[datazonecode].MedianGHIWeek;
                    SDZ.Attributes["income_mean"] = 52* ScotlandIncomes[datazonecode].MeanGHIWeek;
                    Assert.IsTrue(SDZ.Attributes.GetType("income_rtn") == typeof(double));

                    if (IndexedCoast.Intersects(TestGeometry))
                    {
                        SDZ.Attributes["coastal"] = true;
                        //CoastalGeographies.Add(SDZ);
                    }
                    AllGeographies.Add(SDZ);
                }
            }

            Console.WriteLine("Calculating French communes that overlap the 'coastal' area.");
            Dictionary<string, FranceIncome> FranceIncomesClean = new Dictionary<string, FranceIncome>(FranceIncomes.Where(x => x.Value.MED17.HasValue));
            double FrenchAverageIncome = FranceIncomesClean.Average(x => x.Value.MED17.Value);
            foreach (Feature commune in FranceCommuneBoundaries)
            {
                string inseecode = (string)commune.Attributes["insee"];
                commune.Attributes.GetNames().ToList().ForEach(x => commune.Attributes.DeleteAttribute(x));
                commune.Attributes.Add("id", inseecode);
                commune.Attributes.Add("country", "France");
                commune.Attributes.Add("income_rtn", 0.0);
                commune.Attributes.Add("income_mean", 0);
                commune.Attributes.Add("income_mdn", 0);
                commune.Attributes.Add("coastal", false);
                Geometry TestGeometry = commune.Geometry;
                if (TestGeometry.IsValid == false)
                {
                    TestGeometry = new BufferOp(commune.Geometry).GetResultGeometry(0);
                    if (TestGeometry.IsValid == false)
                    {
                        InvalidPolygons.Add(commune);
                    }
                }

                if (FranceIncomesClean.ContainsKey(inseecode))
                {
                    FranceIncome FI = FranceIncomesClean[inseecode];
                    double income = FI.MED17.Value;
                    commune.Attributes["income_rtn"] = (double)(100 * income / FrenchAverageIncome);
                    commune.Attributes["income_mdn"] = FI.MED17;
                    Assert.IsTrue(commune.Attributes.GetType("income_rtn") == typeof(double));

                    if (TestGeometry.IsValid && IndexedCoast.Intersects(TestGeometry))
                    {
                        commune.Attributes["coastal"] = true;
                        //CoastalGeographies.Add(commune);
                    }
                    AllGeographies.Add(commune);
                }
            }
            
            // create the coastal coastal shapefile
            //File.WriteAllText($"CoastalGeographies.geojson", writer.Write(CoastalGeographies));

            // calculate per-country quintiles
            List<string> UniqueCountries = AllGeographies.Select(x => (string)x.Attributes["country"]).Distinct().ToList();
            foreach(string country in UniqueCountries)
            {
                Console.WriteLine($"Calculating income quintiles for {country}.");
                List<double> OrderedListOfAllIncomesInCountry = AllGeographies.Where(x => (string)x.Attributes["country"] == country).Select(x => (double)x.Attributes["income_rtn"]).OrderBy(x => x).ToList();
                Dictionary<double, int> DictionaryOfOrderAllIncomesInCountry = new Dictionary<double, int>();
                for (int i = 0; i < OrderedListOfAllIncomesInCountry.Count; i++)
                {
                    if (DictionaryOfOrderAllIncomesInCountry.ContainsKey(OrderedListOfAllIncomesInCountry[i]) == false)
                    {
                        DictionaryOfOrderAllIncomesInCountry.Add(OrderedListOfAllIncomesInCountry[i], i);
                    }
                }

                List<string> IDsInCountry = AllGeographies.Where(x => (string)x.Attributes["country"] == country).Select(x => (string)x.Attributes["id"]).ToList();
                foreach (string id in IDsInCountry)
                {
                    IFeature elementwithid = AllGeographies.Where(x => (string)x.Attributes["id"] == id).First();
                    double income_rtn = (double)elementwithid.Attributes["income_rtn"];
                    int incomequintile = CalculateNTile(income_rtn, DictionaryOfOrderAllIncomesInCountry, OrderedListOfAllIncomesInCountry.Count, 5);
                    elementwithid.Attributes.Add("incquintile", incomequintile);
                }
            }

            // Create the all geographies shapefile. It is essential that all elements to be output have consistent features.
            using (var streamWriter = new StreamWriter("AllGeographies.prj"))
            {
                streamWriter.Write(GeographicCoordinateSystem.WGS84.WKT);
            }
            GeometryFactory AllGeomFactory = GeometryFactory.Default;
            ShapefileDataWriter Allshpwriter = new ShapefileDataWriter("AllGeographies", AllGeomFactory);
            Allshpwriter.Header = ShapefileDataWriter.GetHeader(AllGeographies[0], AllGeographies.Count);
            Allshpwriter.Write(AllGeographies);


            using (var streamWriter = new StreamWriter("CoastGeographies.prj"))
            {
                streamWriter.Write(GeographicCoordinateSystem.WGS84.WKT);
            }
            FeatureCollection CoastGeographies = new FeatureCollection();
            foreach (Feature fe in AllGeographies)
            {
                if ((bool)fe.Attributes["coastal"] == true)
                {
                    CoastGeographies.Add(fe);
                }
            }
            GeometryFactory CoastGeomFactory = GeometryFactory.Default;
            ShapefileDataWriter Coastshpwriter = new ShapefileDataWriter("CoastGeographies", CoastGeomFactory);
            Coastshpwriter.Header = ShapefileDataWriter.GetHeader(CoastGeographies[0], CoastGeographies.Count);
            Coastshpwriter.Write(CoastGeographies);

            File.WriteAllText($"InvalidBoundaries.geojson", writer.Write(InvalidPolygons));

            // write out a copy of the shapefile attributes table as a CSV
            List<AttributesTableElements> AttributesForCSV = new List<AttributesTableElements>();
            foreach(Feature fe in AllGeographies)
            {
                
                AttributesTableElements ate = new AttributesTableElements();
                ate.coastal = (bool)fe.Attributes["coastal"];
                ate.country = (string)fe.Attributes["country"];
                ate.id = (string)fe.Attributes["id"];
                ate.income_rtn = (int)Math.Round((double)fe.Attributes["income_rtn"],1);
                ate.income_rtnd = (double)fe.Attributes["income_rtn"];
                ate.income_mean = (int)fe.Attributes["income_mean"];
                ate.income_mdn = (int)fe.Attributes["income_mdn"];
                ate.incquintile = (int)fe.Attributes["incquintile"];
                AttributesForCSV.Add(ate);
            }

            // this method of writing the CSV is long -- but it create CSVs better suited for opening in Excel and PowerBI
            if (File.Exists(@"GeographiesAttributesTable.csv"))
            {
                File.Delete(@"GeographiesAttributesTable.csv");
            }
            using (StreamWriter _textWriter = new StreamWriter(File.OpenWrite(@"GeographiesAttributesTable.csv"), new UTF8Encoding(true)))
            {
                using (CsvWriter _csvwriter = new CsvWriter(_textWriter, CultureInfo.CurrentCulture))
                {
                    _csvwriter.WriteRecords(AttributesForCSV);
                }
            }

        }

        public static Geometry Transform(Geometry geom, MathTransform transform)
        {
            geom = geom.Copy();
            geom.Apply(new MTF(transform));
            return geom;
        }

        public static int CalculateNTile(double number, Dictionary<double, int> SortedListOfAllIncomesInCountry, double count, int n)
        {
            // this is an approximation that will only work if allNumbers.Count is much (min. 10x) larger than n!!!
            int position = SortedListOfAllIncomesInCountry[number];
            int ntile = (int)Math.Round(0.5 + ((double)(n * position) / (double)count), 0);
            if (ntile == 0)
            {
                ntile = 1;
            }
            return ntile;
        }
    }

    public class AttributesTableElements
    {
        public string country { get; set; }
        public string id { get; set; }
        public int income_rtn { get; set; }
        public double income_rtnd { get; set; }
        public int income_mean { get; set; }
        public int income_mdn { get; set; }
        public bool coastal { get; set; }
        public int incquintile { get; set; }
    }

    public class ROIIncome
    {
        public string ED { get; set; }
        public int MedianIncome { get; set; }
    }

    public class ScotlandIncome
    {
        public string S2011DataZonecode { get; set; }
        public int MeanGHIWeek { get; set; }
        public int MedianGHIWeek { get; set; }
    }

    public class EandWIncome
    {
        public string MSOAcode { get; set; }
        public double NAIBHC { get; set; }
    }

    public class BelgiumIncome
    {
        public string Code { get; set; }
        public int? Revenumoyenpardec { get; set; }
        public int? Revenumedianpardec { get; set; }
    }

    public class FranceIncome
    {
        public string CODGEO { get; set; }
        public int? NBMENFISC17 { get; set; }
        public int? NBPERSMENFISC17 { get; set; }
        public int? MED17 { get; set; }
        public string PIMP17 { get; set; }
        public string TP6017 { get; set; }
        public string TP60AGE117 { get; set; }
        public string TP60AGE217 { get; set; }
        public string TP60AGE317 { get; set; }
        public string TP60AGE417 { get; set; }
        public string TP60AGE517 { get; set; }
        public string TP60AGE617 { get; set; }
        public string TP60TOL117 { get; set; }
        public string TP60TOL217 { get; set; }
        public string PACT17 { get; set; }
        public string PTSA17 { get; set; }
        public string PCHO17 { get; set; }
        public string PBEN17 { get; set; }
        public string PPEN17 { get; set; }
        public string PPAT17 { get; set; }
        public string PPSOC17 { get; set; }
        public string PPFAM17 { get; set; }
        public string PPMINI17 { get; set; }
        public string PPLOGT17 { get; set; }
        public string PIMPOT17 { get; set; }
        public int? D117 { get; set; }
        public int? D917 { get; set; }
        public string RD17 { get; set; }
    }

    public class NIIncome
    {
        public string SOA2001 { get; set; }
        public string SOA2001_name { get; set; }
        public int IncomeDomainRank_1mostdeprived { get; set; }
    }

    sealed class MTF : NetTopologySuite.Geometries.ICoordinateSequenceFilter
    {
        private readonly MathTransform _mathTransform;

        public MTF(MathTransform mathTransform) => _mathTransform = mathTransform;

        public bool Done => false;
        public bool GeometryChanged => true;
        public void Filter(CoordinateSequence seq, int i)
        {
            double x = seq.GetX(i);
            double y = seq.GetY(i);
            double[] transformed = _mathTransform.Transform(new double[] { x, y });
            seq.SetX(i, transformed[0]);
            seq.SetY(i, transformed[1]);
        }
    }
}