Mirrors/GeoUK
1
0
Fork 0
mirror of https://github.com/IeuanWalker/GeoUK.git synced 2025-10-25 15:19:15 +00:00
Code Issues Packages Projects Releases Wiki Activity

Code refactored

Browse Source
This commit is contained in:
Ieuan Walker
2019-07-12 20:36:23 +01:00
parent 0b8b41a26f
commit c7d1782ab3
24 changed files with 776 additions and 1070 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
GeoUK.OSTN.Tests/ExtensionMethods.cs
View File

@@ -5,15 +5,9 @@ namespace GeoUK.OSTN.Tests
// Source: https://scottlilly.com/c-tip-how-to-check-if-two-double-values-are-equal/
public static class ExtensionMethods
{
public static bool IsApproximatelyEqualTo(this double initialValue, double value)
{
return IsApproximatelyEqualTo(initialValue, value, 0.00001);
}
public static bool IsApproximatelyEqualTo(this double initialValue, double value) => IsApproximatelyEqualTo(initialValue, value, 0.00001);
public static bool IsApproximatelyEqualTo(this double initialValue, double value, double maximumDifferenceAllowed)
{
// Handle comparisons of floating point values that may not be exactly the same
return (Math.Abs(initialValue - value) < maximumDifferenceAllowed);
}
// Handle comparisons of floating point values that may not be exactly the same
public static bool IsApproximatelyEqualTo(this double initialValue, double value, double maximumDifferenceAllowed) => (Math.Abs(initialValue - value) < maximumDifferenceAllowed);
}
}

98
GeoUK.OSTN.Tests/TransformTests.cs
View File

@@ -23,18 +23,17 @@ namespace GeoUK.OSTN.Tests
string line;
while ((line = inputFile.ReadLine()) != null)
{
if (!string.IsNullOrEmpty(line) && line.StartsWith("TP"))
if (string.IsNullOrEmpty(line) || !line.StartsWith("TP")) continue;
string[] values = line.Split(',');
DataPoint point = new DataPoint
{
string[] values = line.Split(',');
DataPoint point = new DataPoint
{
PointID = values[0],
X = double.Parse(values[1]),
Y = double.Parse(values[2]),
Height = double.Parse(values[3])
};
inputData.Add(point);
}
PointID = values[0],
X = double.Parse(values[1]),
Y = double.Parse(values[2]),
Height = double.Parse(values[3])
};
inputData.Add(point);
}
}
@@ -43,24 +42,23 @@ namespace GeoUK.OSTN.Tests
string line;
while ((line = outputFile.ReadLine()) != null)
{
if (!string.IsNullOrEmpty(line) && line.StartsWith("TP"))
if (string.IsNullOrEmpty(line) || !line.StartsWith("TP")) continue;
string[] values = line.Split(',');
DataPoint point = new DataPoint
{
string[] values = line.Split(',');
DataPoint point = new DataPoint
{
PointID = values[0],
X = double.Parse(values[1]),
Y = double.Parse(values[2]),
Height = double.Parse(values[3])
};
outputData[point.PointID] = point;
}
PointID = values[0],
X = double.Parse(values[1]),
Y = double.Parse(values[2]),
Height = double.Parse(values[3])
};
outputData[point.PointID] = point;
}
}
foreach (DataPoint dataPoint in inputData)
{
Osgb36 transformation = Transform.Etrs89ToOsgb(new LatitudeLongitude(dataPoint.X, dataPoint.Y, dataPoint.Height), OstnVersionEnum.OSTN15);
Osgb36 transformation = Transform.Etrs89ToOsgb(new LatitudeLongitude(dataPoint.X, dataPoint.Y, dataPoint.Height));
// Comparing values with a precision of 3 decimals, as they are given in the output file.
bool latitudesEqual = outputData[dataPoint.PointID].X
@@ -90,18 +88,17 @@ namespace GeoUK.OSTN.Tests
string line;
while ((line = inputFile.ReadLine()) != null)
{
if (!string.IsNullOrEmpty(line) && line.StartsWith("TP"))
if (string.IsNullOrEmpty(line) || !line.StartsWith("TP")) continue;
string[] values = line.Split(',');
DataPoint point = new DataPoint
{
string[] values = line.Split(',');
DataPoint point = new DataPoint
{
PointID = values[0],
X = double.Parse(values[1]),
Y = double.Parse(values[2]),
Height = double.Parse(values[3])
};
inputData.Add(point);
}
PointID = values[0],
X = double.Parse(values[1]),
Y = double.Parse(values[2]),
Height = double.Parse(values[3])
};
inputData.Add(point);
}
}
@@ -110,21 +107,19 @@ namespace GeoUK.OSTN.Tests
string line;
while ((line = outputFile.ReadLine()) != null)
{
if (!string.IsNullOrEmpty(line) && line.StartsWith("TP"))
if (string.IsNullOrEmpty(line) || !line.StartsWith("TP")) continue;
string[] values = line.Split(',');
if (values[1] != "RESULT") continue;
DataPoint point = new DataPoint
{
string[] values = line.Split(',');
if (values[1] == "RESULT")
{
DataPoint point = new DataPoint
{
PointID = values[0],
X = double.Parse(values[2]),
Y = double.Parse(values[3]),
Height = double.Parse(values[4])
};
outputData[point.PointID] = point;
}
}
PointID = values[0],
X = double.Parse(values[2]),
Y = double.Parse(values[3]),
Height = double.Parse(values[4])
};
outputData[point.PointID] = point;
}
}
@@ -138,15 +133,8 @@ namespace GeoUK.OSTN.Tests
bool longitudesEqual = outputData[dataPoint.PointID].Y
.IsApproximatelyEqualTo(transformation.Longitude, 0.0000000001);
// Comparing heights with a precision of 4 decimals, as they are given in the output file
// Not implemented
//var heightsEqual = outputData[dataPoint.PointID].Height
// .IsApproximatelyEqualTo(transformation.ElipsoidalHeight, 0.0001);
Assert.True(latitudesEqual);
Assert.True(longitudesEqual);
// Not implemented
//Assert.True(heightsEqual);
}
}
}

2
GeoUK.OSTN.sln.DotSettings Normal file
View File

@@ -0,0 +1,2 @@
<wpf:ResourceDictionary xml:space="preserve" xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:s="clr-namespace:System;assembly=mscorlib" xmlns:ss="urn:shemas-jetbrains-com:settings-storage-xaml" xmlns:wpf="http://schemas.microsoft.com/winfx/2006/xaml/presentation">
<s:Boolean x:Key="/Default/UserDictionary/Words/=northing/@EntryIndexedValue">True</s:Boolean></wpf:ResourceDictionary>

57
GeoUK.OSTN/ResourceManager.cs
View File

@@ -10,27 +10,11 @@ namespace GeoUK.OSTN
{
private static Dictionary<int, OstnDataRecord> _ostn02Data;
public static Dictionary<int, OstnDataRecord> Ostn02Data
{
get
{
if (_ostn02Data == null)
_ostn02Data = RetrieveEmbeddedOSTN(OstnVersionEnum.OSTN02);
return _ostn02Data;
}
}
public static Dictionary<int, OstnDataRecord> Ostn02Data => _ostn02Data ?? (_ostn02Data = RetrieveEmbeddedOSTN(OstnVersionEnum.OSTN02));
private static Dictionary<int, OstnDataRecord> _ostn15Data;
public static Dictionary<int, OstnDataRecord> Ostn15Data
{
get
{
if (_ostn15Data == null)
_ostn15Data = RetrieveEmbeddedOSTN(OstnVersionEnum.OSTN15);
return _ostn15Data;
}
}
public static Dictionary<int, OstnDataRecord> Ostn15Data => _ostn15Data ?? (_ostn15Data = RetrieveEmbeddedOSTN(OstnVersionEnum.OSTN15));
/// <summary>
/// Loads the OSTN data into memory.
@@ -56,11 +40,11 @@ namespace GeoUK.OSTN
switch (ostnVersion)
{
case OstnVersionEnum.OSTN02:
stream = ResourceManager.GetEmbeddedResourceStream(typeof(Transform).GetTypeInfo().Assembly, "OSTN02_OSGM02_GB.txt");
stream = GetEmbeddedResourceStream(typeof(Transform).GetTypeInfo().Assembly, "OSTN02_OSGM02_GB.txt");
break;
case OstnVersionEnum.OSTN15:
stream = ResourceManager.GetEmbeddedResourceStream(typeof(Transform).GetTypeInfo().Assembly, "OSTN15_OSGM15_DataFile.txt");
stream = GetEmbeddedResourceStream(typeof(Transform).GetTypeInfo().Assembly, "OSTN15_OSGM15_DataFile.txt");
break;
default:
@@ -76,21 +60,20 @@ namespace GeoUK.OSTN
string line;
while ((line = reader.ReadLine()) != null)
{
if (!String.IsNullOrWhiteSpace(line))
if (string.IsNullOrWhiteSpace(line)) continue;
string[] values = line.Split(',');
OstnDataRecord record = new OstnDataRecord
{
string[] values = line.Split(',');
OstnDataRecord record = new OstnDataRecord
{
Point_ID = Int32.Parse(values[0]),
ETRS89_Easting = Double.Parse(values[1]),
ETRS89_Northing = Double.Parse(values[2]),
ETRS89_OSGB36_EShift = Double.Parse(values[3]),
ETRS89_OSGB36_NShift = Double.Parse(values[4]),
ETRS89_ODN_HeightShift = Double.Parse(values[5]),
Height_Datum_Flag = Double.Parse(values[6]),
};
data[record.Point_ID] = record;
}
Point_ID = int.Parse(values[0]),
ETRS89_Easting = double.Parse(values[1]),
ETRS89_Northing = double.Parse(values[2]),
ETRS89_OSGB36_EShift = double.Parse(values[3]),
ETRS89_OSGB36_NShift = double.Parse(values[4]),
ETRS89_ODN_HeightShift = double.Parse(values[5]),
Height_Datum_Flag = double.Parse(values[6]),
};
data[record.Point_ID] = record;
}
}
@@ -113,12 +96,12 @@ namespace GeoUK.OSTN
if (!resourcePaths.Any())
{
throw new Exception(String.Format("Resource ending with {0} not found.", resourceFileName));
throw new Exception($"Resource ending with {resourceFileName} not found.");
}
if (resourcePaths.Count() > 1)
if (resourcePaths.Length > 1)
{
throw new Exception(String.Format("Multiple resources ending with {0} found: {1}{2}", resourceFileName, Environment.NewLine, String.Join(Environment.NewLine, resourcePaths)));
throw new Exception($"Multiple resources ending with {resourceFileName} found: {Environment.NewLine}{string.Join(Environment.NewLine, resourcePaths)}");
}
return assembly.GetManifestResourceStream(resourcePaths.Single());

248
GeoUK.OSTN/Transform.cs
View File

@@ -6,157 +6,155 @@ using System.Collections.Generic;
namespace GeoUK.OSTN
{
public static class Transform
{
/// <summary>
/// Loads the OSTN data into memory.
/// </summary>
/// <param name="ostnVersion">If not provided, it will load both OSTN02 and OSTN15 data.</param>
public static void PreloadResources(OstnVersionEnum? ostnVersion = null)
{
ResourceManager.LoadResources(ostnVersion);
}
public static class Transform
{
/// <summary>
/// Loads the OSTN data into memory.
/// </summary>
/// <param name="ostnVersion">If not provided, it will load both OSTN02 and OSTN15 data.</param>
public static void PreloadResources(OstnVersionEnum? ostnVersion = null) => ResourceManager.LoadResources(ostnVersion);
/// <summary>
/// Performs an ETRS89 to OSGB36/ODN datum transformation. Accuracy is approximately 10 centimeters.
/// Whilst very accurate this method is much slower than the Helmert transformation.
/// </summary>
public static Osgb36 Etrs89ToOsgb(LatitudeLongitude coordinates, OstnVersionEnum ostnVersion = OstnVersionEnum.OSTN15)
{
EastingNorthing enCoordinates = Convert.ToEastingNorthing(new Grs80(), new BritishNationalGrid(), coordinates);
return Etrs89ToOsgb(enCoordinates, coordinates.ElipsoidalHeight);
}
/// <summary>
/// Performs an ETRS89 to OSGB36/ODN datum transformation. Accuracy is approximately 10 centimeters.
/// Whilst very accurate this method is much slower than the Helmert transformation.
/// </summary>
public static Osgb36 Etrs89ToOsgb(LatitudeLongitude coordinates)
{
EastingNorthing enCoordinates = Convert.ToEastingNorthing(new Grs80(), new BritishNationalGrid(), coordinates);
return Etrs89ToOsgb(enCoordinates, coordinates.EllipsoidalHeight);
}
private static Osgb36 Etrs89ToOsgb(EastingNorthing coordinates, double ellipsoidHeight, OstnVersionEnum ostnVersion = OstnVersionEnum.OSTN15)
{
Shifts shifts = GetShifts(coordinates, ellipsoidHeight, ostnVersion);
private static Osgb36 Etrs89ToOsgb(EastingNorthing coordinates, double ellipsoidHeight, OstnVersionEnum ostnVersion = OstnVersionEnum.OSTN15)
{
Shifts shifts = GetShifts(coordinates, ellipsoidHeight, ostnVersion);
double easting = coordinates.Easting + shifts.Se;
double northing = coordinates.Northing + shifts.Sn;
double height = ellipsoidHeight - shifts.Sg;
double easting = coordinates.Easting + shifts.Se;
double northing = coordinates.Northing + shifts.Sn;
double height = ellipsoidHeight - shifts.Sg;
return new Osgb36(easting, northing, height, shifts.GeoidDatum);
}
return new Osgb36(easting, northing, height, shifts.GeoidDatum);
}
/// <summary>
/// Performs an OSGB36/ODN to ETRS89 datum transformation. Accuracy is approximately 10 centimeters.
/// Whilst very accurate this method is much slower than the Helmert transformation.
/// </summary>
public static LatitudeLongitude OsgbToEtrs89(Osgb36 coordinates, OstnVersionEnum ostnVersion = OstnVersionEnum.OSTN15)
{
//calculate shifts from OSGB36 point
double errorN = double.MaxValue;
double errorE = double.MaxValue;
EastingNorthing enCoordinates = null;
/// <summary>
/// Performs an OSGB36/ODN to ETRS89 datum transformation. Accuracy is approximately 10 centimeters.
/// Whilst very accurate this method is much slower than the Helmert transformation.
/// </summary>
public static LatitudeLongitude OsgbToEtrs89(Osgb36 coordinates, OstnVersionEnum ostnVersion = OstnVersionEnum.OSTN15)
{
//calculate shifts from OSGB36 point
double errorN = double.MaxValue;
double errorE = double.MaxValue;
EastingNorthing enCoordinates = null;
Shifts shiftsA = GetShifts(coordinates, coordinates.Height, ostnVersion);
Shifts shiftsA = GetShifts(coordinates, coordinates.Height, ostnVersion);
//0.0001 error meters
int iter = 0;
while ((errorN > 0.0001 || errorE > 0.0001) && iter < 10)
{
enCoordinates = new EastingNorthing(coordinates.Easting - shiftsA.Se, coordinates.Northing - shiftsA.Sn);
Shifts shiftsB = GetShifts(enCoordinates, coordinates.Height, ostnVersion);
//0.0001 error meters
int iter = 0;
while ((errorN > 0.0001 || errorE > 0.0001) && iter < 10)
{
enCoordinates = new EastingNorthing(coordinates.Easting - shiftsA.Se, coordinates.Northing - shiftsA.Sn);
Shifts shiftsB = GetShifts(enCoordinates, coordinates.Height, ostnVersion);
errorE = Math.Abs(shiftsA.Se - shiftsB.Se);
errorN = Math.Abs(shiftsA.Sn - shiftsB.Sn);
errorE = Math.Abs(shiftsA.Se - shiftsB.Se);
errorN = Math.Abs(shiftsA.Sn - shiftsB.Sn);
shiftsA = shiftsB;
iter++;
}
shiftsA = shiftsB;
iter++;
}
return Convert.ToLatitudeLongitude(new Wgs84(), new BritishNationalGrid(), enCoordinates);
}
return Convert.ToLatitudeLongitude(new Wgs84(), new BritishNationalGrid(), enCoordinates);
}
private static Shifts GetShifts(EastingNorthing coordinates, double ellipsoidHeight, OstnVersionEnum ostnVersion)
{
//See OS Document: Transformations and OSGM02/OSGM15 user guide chapter 3
Dictionary<int, OstnDataRecord> ostnData = GetOstnData(ostnVersion);
private static Shifts GetShifts(EastingNorthing coordinates, double ellipsoidHeight, OstnVersionEnum ostnVersion)
{
//See OS Document: Transformations and OSGM02/OSGM15 user guide chapter 3
Dictionary<int, OstnDataRecord> ostnData = GetOstnData(ostnVersion);
List<int> recordNumbers = new List<int>();
OstnDataRecord[] records = new OstnDataRecord[4];
List<int> recordNumbers = new List<int>();
OstnDataRecord[] records = new OstnDataRecord[4];
//determine record numbers
int eastIndex = (int)(coordinates.Easting / 1000.0);
int northIndex = (int)(coordinates.Northing / 1000.0);
//determine record numbers
int eastIndex = (int)(coordinates.Easting / 1000.0);
int northIndex = (int)(coordinates.Northing / 1000.0);
double x0 = eastIndex * 1000;
double y0 = northIndex * 1000;
double x0 = eastIndex * 1000;
double y0 = northIndex * 1000;
//work out the four records
recordNumbers.Add(CalculateRecordNumber(eastIndex, northIndex));
recordNumbers.Add(CalculateRecordNumber(eastIndex + 1, northIndex));
recordNumbers.Add(CalculateRecordNumber(eastIndex + 1, northIndex + 1));
recordNumbers.Add(CalculateRecordNumber(eastIndex, northIndex + 1));
//work out the four records
recordNumbers.Add(CalculateRecordNumber(eastIndex, northIndex));
recordNumbers.Add(CalculateRecordNumber(eastIndex + 1, northIndex));
recordNumbers.Add(CalculateRecordNumber(eastIndex + 1, northIndex + 1));
recordNumbers.Add(CalculateRecordNumber(eastIndex, northIndex + 1));
// Get the corresponding reccords from the data dictionary
for (int index = 0; index < 4; index++)
{
records[index] = ostnData[recordNumbers[index]];
}
// Get the corresponding records from the data dictionary
for (int index = 0; index < 4; index++)
{
records[index] = ostnData[recordNumbers[index]];
}
//populate the properties
double se0 = System.Convert.ToDouble(records[0].ETRS89_OSGB36_EShift);
double se1 = System.Convert.ToDouble(records[1].ETRS89_OSGB36_EShift);
double se2 = System.Convert.ToDouble(records[2].ETRS89_OSGB36_EShift);
double se3 = System.Convert.ToDouble(records[3].ETRS89_OSGB36_EShift);
//populate the properties
double se0 = System.Convert.ToDouble(records[0].ETRS89_OSGB36_EShift);
double se1 = System.Convert.ToDouble(records[1].ETRS89_OSGB36_EShift);
double se2 = System.Convert.ToDouble(records[2].ETRS89_OSGB36_EShift);
double se3 = System.Convert.ToDouble(records[3].ETRS89_OSGB36_EShift);
double sn0 = System.Convert.ToDouble(records[0].ETRS89_OSGB36_NShift);
double sn1 = System.Convert.ToDouble(records[1].ETRS89_OSGB36_NShift);
double sn2 = System.Convert.ToDouble(records[2].ETRS89_OSGB36_NShift);
double sn3 = System.Convert.ToDouble(records[3].ETRS89_OSGB36_NShift);
double sn0 = System.Convert.ToDouble(records[0].ETRS89_OSGB36_NShift);
double sn1 = System.Convert.ToDouble(records[1].ETRS89_OSGB36_NShift);
double sn2 = System.Convert.ToDouble(records[2].ETRS89_OSGB36_NShift);
double sn3 = System.Convert.ToDouble(records[3].ETRS89_OSGB36_NShift);
double sg0 = System.Convert.ToDouble(records[0].ETRS89_ODN_HeightShift);
double sg1 = System.Convert.ToDouble(records[1].ETRS89_ODN_HeightShift);
double sg2 = System.Convert.ToDouble(records[2].ETRS89_ODN_HeightShift);
double sg3 = System.Convert.ToDouble(records[3].ETRS89_ODN_HeightShift);
double sg0 = System.Convert.ToDouble(records[0].ETRS89_ODN_HeightShift);
double sg1 = System.Convert.ToDouble(records[1].ETRS89_ODN_HeightShift);
double sg2 = System.Convert.ToDouble(records[2].ETRS89_ODN_HeightShift);
double sg3 = System.Convert.ToDouble(records[3].ETRS89_ODN_HeightShift);
double dx = coordinates.Easting - x0;
double dy = coordinates.Northing - y0;
double dx = coordinates.Easting - x0;
double dy = coordinates.Northing - y0;
double t = dx / 1000.0;
double u = dy / 1000.0;
double t = dx / 1000.0;
double u = dy / 1000.0;
Shifts shifts = new Shifts();
Shifts shifts = new Shifts
{
Se = (1 - t) * (1 - u) * se0 + t * (1 - u) * se1 + t * u * se2 + (1 - t) * u * se3,
Sn = (1 - t) * (1 - u) * sn0 + t * (1 - u) * sn1 + t * u * sn2 + (1 - t) * u * sn3,
Sg = (1 - t) * (1 - u) * sg0 + t * (1 - u) * sg1 + t * u * sg2 + (1 - t) * u * sg3,
shifts.Se = (1 - t) * (1 - u) * se0 + t * (1 - u) * se1 + t * u * se2 + (1 - t) * u * se3;
shifts.Sn = (1 - t) * (1 - u) * sn0 + t * (1 - u) * sn1 + t * u * sn2 + (1 - t) * u * sn3;
shifts.Sg = (1 - t) * (1 - u) * sg0 + t * (1 - u) * sg1 + t * u * sg2 + (1 - t) * u * sg3;
GeoidDatum = (Osgb36GeoidDatum)System.Convert.ToInt32(records[0].Height_Datum_Flag)
};
shifts.GeoidDatum = (Osgb36GeoidDatum)System.Convert.ToInt32(records[0].Height_Datum_Flag);
return shifts;
}
return shifts;
}
/// <summary>
/// Calculates a data file record number.
/// </summary>
/// <param name="eastIndex"></param>
/// <param name="northIndex"></param>
/// <returns></returns>
private static int CalculateRecordNumber(int eastIndex, int northIndex)
{
return eastIndex + (northIndex * 701) + 1;
}
/// <summary>
/// Calculates a data file record number.
/// </summary>
/// <param name="eastIndex"></param>
/// <param name="northIndex"></param>
/// <returns></returns>
private static int CalculateRecordNumber(int eastIndex, int northIndex)
{
return eastIndex + (northIndex * 701) + 1;
}
/// <summary>
/// Retrieves the parsed OSTN data from the embedded resource file.
/// </summary>
/// <param name="ostnVersion"></param>
/// <returns></returns>
private static Dictionary<int, OstnDataRecord> GetOstnData(OstnVersionEnum ostnVersion)
{
switch (ostnVersion)
{
case OstnVersionEnum.OSTN02:
return ResourceManager.Ostn02Data;
/// <summary>
/// Retrieves the parsed OSTN data from the embedded resource file.
/// </summary>
/// <param name="ostnVersion"></param>
/// <returns></returns>
private static Dictionary<int, OstnDataRecord> GetOstnData(OstnVersionEnum ostnVersion)
{
switch (ostnVersion)
{
case OstnVersionEnum.OSTN02:
return ResourceManager.Ostn02Data;
case OstnVersionEnum.OSTN15:
return ResourceManager.Ostn15Data;
case OstnVersionEnum.OSTN15:
return ResourceManager.Ostn15Data;
default:
throw new NotImplementedException();
}
}
}
default:
throw new NotImplementedException();
}
}
}
}

156
GeoUK/Convert.cs
View File

@@ -19,11 +19,11 @@ namespace GeoUK
/// <param name = "ellipsoid"></param>
/// <param name="projection">Projection.</param>
/// <param name="coordinates">Coordinates.</param>
public static LatitudeLongitude ToLatitudeLongitude(Ellipsoid ellipsoid, Projections.Projection projection, EastingNorthing coordinates)
public static LatitudeLongitude ToLatitudeLongitude(Ellipsoid ellipsoid, Projection projection, EastingNorthing coordinates)
{
double M;
double N = coordinates.Northing;
double E = coordinates.Easting;
double N = coordinates.Northing;
double E = coordinates.Easting;
//from OS Guide
//constants needed are a Semi-Major Axis , b , e2 , N0 , E0 , F0 ,φ0 , and λ0
@@ -35,25 +35,8 @@ namespace GeoUK
double lon0 = ToRadians(projection.TrueOriginLongitude);
double E0 = projection.TrueOriginEasting;
double N0 = projection.TrueOriginNorthing;
double lat = ((N - N0) / (a * F0)) + lat0;
/*
Console.WriteLine ("a: {0}", a);
Console.WriteLine ("b: {0}", b);
Console.WriteLine ("e2: {0}", e2);
Console.WriteLine ("F0: {0}", F0);
Console.WriteLine ("lat0 (deg): {0}", projection.TrueOriginLatitude);
Console.WriteLine ("lon0 (deg): {0}", projection.TrueOriginLongitude);
Console.WriteLine ("lat0 (rad): {0}", lat0);
Console.WriteLine ("lon0 (rad): {0}", lon0);
Console.WriteLine ("E0: {0}", E0);
Console.WriteLine ("N0: {0}", N0);
Console.WriteLine ("N: {0}", N);
Console.WriteLine ("E: {0}", E);
*/
double lat = ((N - N0) / (a * F0)) + lat0;
//double n = (a - b) / (a + b);
//check for error and reiterate as required
int loopCount = 0;
do
@@ -62,12 +45,6 @@ namespace GeoUK
lat += ((N - N0 - M) / (a * F0));
/*
Console.WriteLine ("\nlat #{0}: {1}", loopCount + 1, lat);
Console.WriteLine ("M #{0}: {1}", loopCount + 1, M);
Console.WriteLine ("N - N0 - M #{0}: {1}\n", loopCount + 1, N - N0 - M);
*/
loopCount++;
} while (!IsNearlyZero(N - N0 - M, 1e-16) & loopCount < 10);
@@ -75,47 +52,32 @@ namespace GeoUK
double p = a * F0 * (1 - e2) * Math.Pow((1 - e2 * Math.Pow(Math.Sin(lat), 2)), -1.5);
double n2 = (v / p) - 1;
double VII = Math.Tan(lat) / (2 * p * v);
double VII = Math.Tan(lat) / (2 * p * v);
double VIIIa = Math.Tan(lat) / (24 * p * Math.Pow(v, 3));
double VIIIb = 5 + (3 * Math.Pow(Math.Tan(lat), 2)) + n2 - 9 * (Math.Pow(Math.Tan(lat), 2)) * n2;
double VIII = VIIIa * VIIIb;
double VIIIa = Math.Tan(lat) / (24 * p * Math.Pow(v, 3));
double VIIIb = 5 + (3 * Math.Pow(Math.Tan(lat), 2)) + n2 - 9 * (Math.Pow(Math.Tan(lat), 2)) * n2;
double VIII = VIIIa * VIIIb;
double IXa = Math.Tan(lat) / (720 * p * Math.Pow(v, 5));
double IXb = 61 + (90 * Math.Pow(Math.Tan(lat), 2)) + (45 * Math.Pow(Math.Tan(lat), 4));
double IX = IXa * IXb;
double IXa = Math.Tan(lat) / (720 * p * Math.Pow(v, 5));
double IXb = 61 + (90 * Math.Pow(Math.Tan(lat), 2)) + (45 * Math.Pow(Math.Tan(lat), 4));
double IX = IXa * IXb;
double X = MathEx.Secant(lat) / v;
double X = MathEx.Secant(lat) / v;
double XIa = MathEx.Secant(lat) / (6 * Math.Pow(v, 3));
double XIb = v / p + (2 * Math.Pow(Math.Tan(lat), 2));
double XI = XIa * XIb;
double XIa = MathEx.Secant(lat) / (6 * Math.Pow(v, 3));
double XIb = v / p + (2 * Math.Pow(Math.Tan(lat), 2));
double XI = XIa * XIb;
double XIIa = MathEx.Secant(lat) / (120 * Math.Pow(v, 5));
double XIIb = 5 + (28 * Math.Pow(Math.Tan(lat), 2)) + (24 * Math.Pow(Math.Tan(lat), 4));
double XII = XIIa * XIIb;
double XIIa = MathEx.Secant(lat) / (120 * Math.Pow(v, 5));
double XIIb = 5 + (28 * Math.Pow(Math.Tan(lat), 2)) + (24 * Math.Pow(Math.Tan(lat), 4));
double XII = XIIa * XIIb;
double XIIAa = MathEx.Secant(lat) / (5040 * Math.Pow(v, 7));
double XIIAb = 61 + (662 * Math.Pow(Math.Tan(lat), 2)) + (1320 * Math.Pow(Math.Tan(lat), 4)) + (720 * Math.Pow(Math.Tan(lat), 6));
double XIIA = XIIAa * XIIAb;
double XIIAa = MathEx.Secant(lat) / (5040 * Math.Pow(v, 7));
double XIIAb = 61 + (662 * Math.Pow(Math.Tan(lat), 2)) + (1320 * Math.Pow(Math.Tan(lat), 4)) + (720 * Math.Pow(Math.Tan(lat), 6));
double XIIA = XIIAa * XIIAb;
lat = lat - VII * Math.Pow(E - E0, 2) + VIII * Math.Pow(E - E0, 4) - IX * Math.Pow(E - E0, 6);
double lon = lon0 + X * (E - E0) - XI * Math.Pow(E - E0, 3) + XII * Math.Pow(E - E0, 5) - XIIA * Math.Pow(E - E0, 7);
/*
Console.WriteLine ("v: {0}", v);
Console.WriteLine ("p: {0}", p);
Console.WriteLine ("n2: {0}", n2);
Console.WriteLine ("VII: {0}", VII);
Console.WriteLine ("VIII: {0}", VIII);
Console.WriteLine ("IX: {0}", IX);
Console.WriteLine ("X: {0}", X);
Console.WriteLine ("XI: {0}", XI);
Console.WriteLine ("XII: {0}", XII);
Console.WriteLine ("XIIA: {0}", XIIA);
Console.WriteLine ("lat (final): {0}", lat);
Console.WriteLine ("lon: {0}", lon);
*/
double lon = lon0 + X * (E - E0) - XI * Math.Pow(E - E0, 3) + XII * Math.Pow(E - E0, 5) - XIIA * Math.Pow(E - E0, 7);
return new LatitudeLongitude(ToDegrees(lat), ToDegrees(lon));
}
@@ -124,7 +86,7 @@ namespace GeoUK
{
if (x < 0)
{
x = x * -1;
x *= -1;
}
return (x < 0 && x > tolerance * -1) || (x >= 0 && x < tolerance);
@@ -148,20 +110,14 @@ namespace GeoUK
/// </summary>
/// <param name="degrees"></param>
/// <returns></returns>
public static double ToRadians(double degrees)
{
return degrees * (Math.PI / 180);
}
public static double ToRadians(double degrees) => degrees * (Math.PI / 180);
/// <summary>
/// Converts radians to decimal degrees.
/// </summary>
/// <param name="radians"></param>
/// <returns></returns>
public static double ToDegrees(double radians)
{
return radians * (180 / Math.PI);
}
public static double ToDegrees(double radians) => radians * (180 / Math.PI);
/// <summary>
/// Converts cartesian coordinates to grid eastings and northings for any Transverse Mercator map projection, including the Ordnance Survey National Grid.
@@ -173,7 +129,7 @@ namespace GeoUK
/// (latitude, longitude), use the GRS80 ellipsoid. Use the same National Grid projection
/// constants for both ETRS89 and OSGB36 coordinates.
/// </remarks>
public static EastingNorthing ToEastingNorthing(Ellipsoid ellipsoid, Projections.Projection projection, Cartesian coordinates)
public static EastingNorthing ToEastingNorthing(Ellipsoid ellipsoid, Projection projection, Cartesian coordinates)
{
LatitudeLongitude coords = ToLatitudeLongitude(ellipsoid, coordinates);
return ToEastingNorthing(ellipsoid, projection, coords);
@@ -193,7 +149,7 @@ namespace GeoUK
/// <param name="projection"></param>
/// <param name="coordinates"></param>
/// <returns></returns>
public static EastingNorthing ToEastingNorthing(Ellipsoid ellipsoid, Projections.Projection projection, LatitudeLongitude coordinates)
public static EastingNorthing ToEastingNorthing(Ellipsoid ellipsoid, Projection projection, LatitudeLongitude coordinates)
{
double lat = ToRadians(coordinates.Latitude);
double lon = ToRadians(coordinates.Longitude);
@@ -221,14 +177,8 @@ namespace GeoUK
//B5
double n2 = v / p - 1;
//B6
double M = CalculateM(lat, lat0, a, b, F0);
// double Ma = (1 + n + ((5.0 / 4.0) * Math.Pow (n, 2)) + ((5.0 / 4.0) * Math.Pow (n, 3))) * (lat - lat0);
// double Mb = ((3 * n) + (3 * Math.Pow (n, 2)) + ((21.0 / 8.0) * Math.Pow (n, 3))) * Math.Sin (lat - lat0) * Math.Cos (lat + lat0);
// double Mc = (((15.0 / 8.0) * Math.Pow (n, 2)) + (15.0 / 8.0) * Math.Pow (n, 3)) * Math.Sin (2 * (lat - lat0)) * Math.Cos (2 * (lat + lat0));
// double Md = (35.0 / 24.0 * Math.Pow (n, 3)) * Math.Sin (3 * (lat - lat0)) * Math.Cos (3 * (lat + lat0));
// double M = b * F0 * (Ma - Mb + Mc - Md);
//B6
double M = CalculateM(lat, lat0, a, b, F0);
double I = M + N0;
double II = (v / 2) * Math.Sin(lat) * Math.Cos(lat);
@@ -244,12 +194,12 @@ namespace GeoUK
//B8
double E = E0 + (IV * (lon - lon0)) + (V * Math.Pow((lon - lon0), 3)) + (VI * Math.Pow((lon - lon0), 5));
return new EastingNorthing(E, N, coordinates.ElipsoidalHeight); //height is still with respect to the ellipsoid
return new EastingNorthing(E, N, coordinates.EllipsoidalHeight); //height is still with respect to the ellipsoid
}
public static Cartesian ToCartesian(Ellipsoid ellipsoid, Projection projection, EastingNorthing coordinates)
{
LatitudeLongitude latLongCoordinates = ToLatitudeLongitude(
LatitudeLongitude latLongCoordinates = ToLatitudeLongitude(
ellipsoid,
projection,
coordinates);
@@ -258,7 +208,7 @@ namespace GeoUK
}
/// <summary>
/// Converts latitude, longitude and elipsoidal height coordinates to cartesian coordinates using the same ellipsoid.
/// Converts latitude, longitude and ellipsoidal height coordinates to cartesian coordinates using the same ellipsoid.
/// Please note this is not a transformation between ellipsoids.
/// </summary>
/// <param name="ellipsoid"></param>
@@ -268,7 +218,7 @@ namespace GeoUK
{
double lat = ToRadians(coordinates.Latitude);
double lon = ToRadians(coordinates.Longitude);
double height = coordinates.ElipsoidalHeight;
double height = coordinates.EllipsoidalHeight;
double e2 = ellipsoid.EccentricitySquared;
double a = ellipsoid.SemiMajorAxis;
@@ -283,7 +233,7 @@ namespace GeoUK
}
/// <summary>
/// Converts cartesian coordinates to latitude, longitude and elipsoidal height using the same ellipsoid.
/// Converts cartesian coordinates to latitude, longitude and ellipsoidal height using the same ellipsoid.
/// Please note this is not a transformation between ellipsoids.
/// </summary>
/// <param name="ellipsoid"></param>
@@ -318,49 +268,38 @@ namespace GeoUK
/// <param name="minutes"></param>
/// <param name="seconds"></param>
/// <returns></returns>
public static double ToDecimelDegrees(int degrees, int minutes, double seconds)
public static double ToDecimalDegrees(int degrees, int minutes, double seconds)
{
//determine seconds as minutes
double m = minutes + (seconds / 60.0);
return ToDecimelDegrees(degrees, m);
return ToDecimalDegrees(degrees, m);
}
/// <summary>
/// Converts degrees and decimel minutes to decimal degrees.
/// Converts degrees and decimal minutes to decimal degrees.
/// </summary>
/// <param name="degrees"></param>
/// <param name="minutes"></param>
/// <returns></returns>
public static double ToDecimelDegrees(int degrees, double minutes)
{
//determine minutes as derees
return degrees + (minutes / 60.0);
}
public static double ToDecimalDegrees(int degrees, double minutes) => degrees + (minutes / 60.0);
private static double Div(double value, double divisor)
{
double dblResult = 0;
bool blnNegative = false;
//make the division
dblResult = value / divisor;
double dblResult = value / divisor;
//do all calculations on positive numbers
bool blnNegative = false;
if (dblResult < 0)
{
blnNegative = true;
dblResult = dblResult * -1;
dblResult *= -1;
}
//see if there is any remainder
if (dblResult % 1 > 0)
{
dblResult = Math.Ceiling(dblResult) - 1;
}
else
{
dblResult = System.Convert.ToInt32(dblResult, CultureInfo.InvariantCulture);
}
dblResult = dblResult % 1 > 0
? Math.Ceiling(dblResult) - 1
: System.Convert.ToInt32(dblResult, CultureInfo.InvariantCulture);
if (blnNegative)
{
@@ -371,13 +310,10 @@ namespace GeoUK
}
/// <summary>
/// Helper funtion to reverse the sign of a value. Helps code to be more readable.
/// Helper function to reverse the sign of a value. Helps code to be more readable.
/// </summary>
/// <param name="value"></param>
/// <returns></returns>
private static double Negate(double value)
{
return value * -1.0;
}
private static double Negate(double value) => value * -1.0;
}
}

24
GeoUK/Coordinates/Cartesian.cs
View File

@@ -5,9 +5,6 @@ namespace GeoUK.Coordinates
/// </summary>
public class Cartesian
{
private double _x;
private double _z;
private double _y;
/// <summary>
/// Constructor.
@@ -17,33 +14,24 @@ namespace GeoUK.Coordinates
/// <param name="z"></param>
public Cartesian(double x, double y, double z)
{
_x = x;
_y = y;
_z = z;
X = x;
Y = y;
Z = z;
}
/// <summary>
/// Returns the X axis parameter.
/// </summary>
public double X
{
get { return _x; }
}
public double X { get; }
/// <summary>
/// Returns the Y axis parameter.
/// </summary>
public double Y
{
get { return _y; }
}
public double Y { get; }
/// <summary>
/// Returns the Z axis parameter.
/// </summary>
public double Z
{
get { return _z; }
}
public double Z { get; }
}
}

32
GeoUK/Coordinates/EastingNorthing.cs
View File

@@ -6,9 +6,6 @@ namespace GeoUK.Coordinates
/// </summary>
public class EastingNorthing
{
private double _easting;
private double _northing;
private double _height;
/// <summary>
/// Constructor.
@@ -17,9 +14,9 @@ namespace GeoUK.Coordinates
/// <param name="northing"></param>
public EastingNorthing(double easting, double northing)
{
_easting = easting;
_northing = northing;
_height = 0;
Easting = easting;
Northing = northing;
Height = 0;
}
/// <summary>
@@ -30,33 +27,24 @@ namespace GeoUK.Coordinates
/// <param name="height"></param>
public EastingNorthing(double easting, double northing, double height)
{
_easting = easting;
_northing = northing;
_height = height;
Easting = easting;
Northing = northing;
Height = height;
}
/// <summary>
/// Retruns the easting parameter.
/// Returns the easting parameter.
/// </summary>
public double Easting
{
get { return _easting; }
}
public double Easting { get; }
/// <summary>
/// returns the northing parameter.
/// </summary>
public double Northing
{
get { return _northing; }
}
public double Northing { get; }
/// <summary>
/// Returns the height parameter.
/// </summary>
public double Height
{
get { return _height; }
}
public double Height { get; }
}
}

103
GeoUK/Coordinates/LatitudeLongitude.cs
View File

@@ -1,71 +1,50 @@
namespace GeoUK.Coordinates
{
/// <summary>
/// This immutable class represents a set of latitude/longitude/ellipsoidal height coordinates.
/// </summary>
public class LatitudeLongitude
/// <summary>
/// This immutable class represents a set of latitude/longitude/ellipsoidal height coordinates.
/// </summary>
public class LatitudeLongitude
{
private double _degreesLatitude = 0.0;
private double _degreesLongitude = 0.0;
private double _elipsoidalHeight = 0.0;
{
/// <summary>
/// Constructor.
/// </summary>
/// <param name="degreesLatitude"></param>
/// <param name="degreesLongitude"></param>
public LatitudeLongitude(double degreesLatitude, double degreesLongitude)
{
_degreesLatitude = degreesLatitude;
_degreesLongitude = degreesLongitude;
_elipsoidalHeight = 0.0;
}
/// <summary>
/// Constructor.
/// </summary>
/// <param name="degreesLatitude"></param>
/// <param name="degreesLongitude"></param>
public LatitudeLongitude(double degreesLatitude, double degreesLongitude)
{
Latitude = degreesLatitude;
Longitude = degreesLongitude;
EllipsoidalHeight = 0.0;
}
/// <summary>
/// Constructor.
/// </summary>
/// <param name="degreesLatitude"></param>
/// <param name="degreesLongitude"></param>
/// <param name="elipsoidalHeight"></param>
public LatitudeLongitude(double degreesLatitude, double degreesLongitude, double elipsoidalHeight)
{
_degreesLatitude = degreesLatitude;
_degreesLongitude = degreesLongitude;
_elipsoidalHeight = elipsoidalHeight;
}
/// <summary>
/// Constructor.
/// </summary>
/// <param name="degreesLatitude"></param>
/// <param name="degreesLongitude"></param>
/// <param name="ellipsoidalHeight"></param>
public LatitudeLongitude(double degreesLatitude, double degreesLongitude, double ellipsoidalHeight)
{
Latitude = degreesLatitude;
Longitude = degreesLongitude;
EllipsoidalHeight = ellipsoidalHeight;
}
/// <summary>
/// Returns latitude in degrees.
/// </summary>
public double Latitude
{
get
{
return _degreesLatitude;
}
}
/// <summary>
/// Returns latitude in degrees.
/// </summary>
public double Latitude { get; }
/// <summary>
/// Returns longitude in degrees.
/// </summary>
public double Longitude
{
get
{
return _degreesLongitude;
}
}
/// <summary>
/// Returns longitude in degrees.
/// </summary>
public double Longitude { get; }
/// <summary>
/// returns elipsoidal height in meters.
/// </summary>
public double ElipsoidalHeight
{
get
{
return _elipsoidalHeight;
}
}
}
/// <summary>
/// returns ellipsoidal height in meters.
/// </summary>
public double EllipsoidalHeight { get; }
}
}

224
GeoUK/Coordinates/Osgb36.cs
View File

@@ -2,140 +2,136 @@ using System;
namespace GeoUK.Coordinates
{
/// <summary>
/// This immutable class, derived from EastingNorthingCoordinates, provides a convenient means
/// to represent OSGB36 eastings and northings.
/// </summary>
/// <remarks>
/// Eastings and northings are represented in British National Grid and Height is specifgied
/// in meters based on the geoid datum returned by the RegionGeoidDatum property.
/// </remarks>
public class Osgb36 : EastingNorthing
{
private Osgb36GeoidDatum _datum = Osgb36GeoidDatum.OutsideModelBoundary;
/// <summary>
/// This immutable class, derived from EastingNorthingCoordinates, provides a convenient means
/// to represent OSGB36 eastings and northings.
/// </summary>
/// <remarks>
/// Eastings and northings are represented in British National Grid and Height is specified
/// in meters based on the geoid datum returned by the RegionGeoidDatum property.
/// </remarks>
public class Osgb36 : EastingNorthing
{
/// <summary>
/// Initializes a new instance of the <see cref="GeoUK.Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="easting">Easting.</param>
/// <param name="northing">Northing.</param>
public Osgb36(double easting, double northing)
: base(easting, northing, 0)
{
_datum = Osgb36GeoidDatum.NewlynUkMainland;
}
/// <summary>
/// Initializes a new instance of the <see cref="Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="easting">Easting.</param>
/// <param name="northing">Northing.</param>
public Osgb36(double easting, double northing)
: base(easting, northing, 0)
{
RegionGeoidDatum = Osgb36GeoidDatum.NewlynUkMainland;
}
/// <summary>
/// Initializes a new instance of the <see cref="GeoUK.Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="eastingNorthingCoordinates">Easting northing coordinates.</param>
public Osgb36(EastingNorthing eastingNorthingCoordinates)
: base(eastingNorthingCoordinates.Easting, eastingNorthingCoordinates.Northing, eastingNorthingCoordinates.Height)
{
_datum = Osgb36GeoidDatum.NewlynUkMainland;
}
/// <summary>
/// Initializes a new instance of the <see cref="Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="eastingNorthingCoordinates">Easting northing coordinates.</param>
public Osgb36(EastingNorthing eastingNorthingCoordinates)
: base(eastingNorthingCoordinates.Easting, eastingNorthingCoordinates.Northing, eastingNorthingCoordinates.Height)
{
RegionGeoidDatum = Osgb36GeoidDatum.NewlynUkMainland;
}
/// <summary>
/// Initializes a new instance of the <see cref="GeoUK.Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="eastingNorthingCoordinates">Easting northing coordinates.</param>
/// <param name="datum">Datum.</param>
public Osgb36(EastingNorthing eastingNorthingCoordinates, Osgb36GeoidDatum datum)
: base(eastingNorthingCoordinates.Easting, eastingNorthingCoordinates.Northing, eastingNorthingCoordinates.Height)
{
_datum = datum;
}
/// <summary>
/// Initializes a new instance of the <see cref="Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="eastingNorthingCoordinates">Easting northing coordinates.</param>
/// <param name="datum">Datum.</param>
public Osgb36(EastingNorthing eastingNorthingCoordinates, Osgb36GeoidDatum datum)
: base(eastingNorthingCoordinates.Easting, eastingNorthingCoordinates.Northing, eastingNorthingCoordinates.Height)
{
RegionGeoidDatum = datum;
}
/// <summary>
/// Initializes a new instance of the <see cref="GeoUK.Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="easting">Easting.</param>
/// <param name="northing">Northing.</param>
/// <param name="height">Height.</param>
/// <param name="datum">Datum.</param>
public Osgb36(double easting, double northing, double height, Osgb36GeoidDatum datum)
: base(easting, northing, height)
{
_datum = datum;
}
/// <summary>
/// Initializes a new instance of the <see cref="Coordinates.Osgb36Cordinates"/> class.
/// </summary>
/// <param name="easting">Easting.</param>
/// <param name="northing">Northing.</param>
/// <param name="height">Height.</param>
/// <param name="datum">Datum.</param>
public Osgb36(double easting, double northing, double height, Osgb36GeoidDatum datum)
: base(easting, northing, height)
{
RegionGeoidDatum = datum;
}
/// <summary>
/// Returns the Local Geoid datum in use. other property values should be
/// considered invalid if this property is set to OutsideModelBoundary.
/// </summary>
public Osgb36GeoidDatum RegionGeoidDatum
{
get { return _datum; }
}
/// <summary>
/// Returns the Local Geoid datum in use. other property values should be
/// considered invalid if this property is set to OutsideModelBoundary.
/// </summary>
public Osgb36GeoidDatum RegionGeoidDatum { get; }
public string MapReference
{
get
{
/*
public string MapReference
{
get
{
/*
10km (2-figure) Grid Reference: SO84 = 380000 Easting 240000 Northing
1km (4-figure) Grid Reference: NS2468 = 224000 Easting 668000 Northing
100m (6-figure) Grid Reference: TL123456 = 512300 Easting 245600 Northing
*/
double easting = this.Easting;
double northing = this.Northing;
*/
double easting = Easting;
double northing = Northing;
string bngSquare = GetBngSquare(easting, northing);
string bngSquare = GetBngSquare(easting, northing);
//get the number of complete 500k squares
int indexNorthing = (int)Math.Floor(northing / 500000);
int indexEasting = (int)Math.Floor(easting / 500000);
//get the number of complete 500k squares
int indexNorthing = (int)Math.Floor(northing / 500000);
int indexEasting = (int)Math.Floor(easting / 500000);
//reduce E and N by the number of 500k squares
northing = northing - indexNorthing * 500000;
easting = easting - indexEasting * 500000;
//reduce E and N by the number of 500k squares
northing -= indexNorthing * 500000;
easting -= indexEasting * 500000;
//reduce by the number of 100k squares within the 500k square.
indexNorthing = (int)Math.Floor(northing) / 100000;
indexEasting = (int)Math.Floor(easting) / 100000;
//reduce by the number of 100k squares within the 500k square.
indexNorthing = (int)Math.Floor(northing) / 100000;
indexEasting = (int)Math.Floor(easting) / 100000;
northing = northing - indexNorthing * 100000;
easting = easting - indexEasting * 100000;
northing -= indexNorthing * 100000;
easting -= indexEasting * 100000;
northing = Math.Round(northing / 100);
easting = Math.Round(easting / 100);
return string.Format("{0}{1}{2}", bngSquare, Math.Round(easting).ToString("000"), Math.Round(northing).ToString("000"));
}
}
northing = Math.Round(northing / 100);
easting = Math.Round(easting / 100);
return $"{bngSquare}{Math.Round(easting):000}{Math.Round(northing):000}";
}
}
/// <summary>
/// Returns the two letter OS code based on easting and northing in metres.
/// </summary>
/// <returns>The square with northing.</returns>
/// <param name="northing">Northing.</param>
/// <param name="easting">Easting.</param>
public static string GetBngSquare(double easting, double northing)
{
string result = string.Empty;
/// <summary>
/// Returns the two letter OS code based on easting and northing in metres.
/// </summary>
/// <returns>The square with northing.</returns>
/// <param name="northing">Northing.</param>
/// <param name="easting">Easting.</param>
public static string GetBngSquare(double easting, double northing)
{
string result = string.Empty;
//test for our upper and lower limits
if (easting >= 0 && easting < 700000 && northing >= 0 && northing < 1300000)
{
char[] firstChar = new char[6] { 'S', 'N', 'H', 'T', 'O', 'J' };
char[] secondChar = new char[25] { 'V', 'Q', 'L', 'F', 'A', 'W', 'R', 'M', 'G', 'B', 'X', 'S', 'N', 'H', 'C', 'Y', 'T', 'O', 'J', 'D', 'Z', 'U', 'P', 'K', 'E' };
//test for our upper and lower limits
if (easting >= 0 && easting < 700000 && northing >= 0 && northing < 1300000)
{
char[] firstChar = new char[6] { 'S', 'N', 'H', 'T', 'O', 'J' };
char[] secondChar = new char[25] { 'V', 'Q', 'L', 'F', 'A', 'W', 'R', 'M', 'G', 'B', 'X', 'S', 'N', 'H', 'C', 'Y', 'T', 'O', 'J', 'D', 'Z', 'U', 'P', 'K', 'E' };
//calculate the first letter
int indexNorthing = (int)Math.Floor(northing / 500000);
int indexEasting = (int)Math.Floor(easting / 500000);
//calculate the first letter
int indexNorthing = (int)Math.Floor(northing / 500000);
int indexEasting = (int)Math.Floor(easting / 500000);
//get the first char
char chr1 = firstChar[(indexEasting * 3) + indexNorthing];
//get the first char
char chr1 = firstChar[(indexEasting * 3) + indexNorthing];
//to get the second letter we subtract the number of 500km sectors calculated above
indexNorthing = (int)Math.Floor((northing - (indexNorthing * 500000)) / 100000);
indexEasting = (int)Math.Floor((easting - (indexEasting * 500000)) / 100000);
//to get the second letter we subtract the number of 500km sectors calculated above
indexNorthing = (int)Math.Floor((northing - (indexNorthing * 500000)) / 100000);
indexEasting = (int)Math.Floor((easting - (indexEasting * 500000)) / 100000);
//get the second char
char chr2 = secondChar[(indexEasting * 5) + indexNorthing];
//get the second char
char chr2 = secondChar[(indexEasting * 5) + indexNorthing];
result = string.Format("{0}{1}", chr1, chr2);
}
return result;
}
}
result = $"{chr1}{chr2}";
}
return result;
}
}
}

30
GeoUK/Ellipsoids/Airy1830.cs
View File

@@ -1,19 +1,19 @@
namespace GeoUK.Ellipsoids
{
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an Airy1930 ellipsoid and is provided for convienience.
/// </summary>
public class Airy1830 : Ellipsoid
{
private const double C_SEMI_MAJOR_AXIS = 6377563.396; //a
private const double C_SEMI_MINOR_AXIS = 6356256.909; //b
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an Airy1930 ellipsoid and is provided for convenience.
/// </summary>
public class Airy1830 : Ellipsoid
{
private const double C_SEMI_MAJOR_AXIS = 6377563.396;
private const double C_SEMI_MINOR_AXIS = 6356256.909;
/// <summary>
/// Constructor.
/// </summary>
public Airy1830()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
/// <summary>
/// Constructor.
/// </summary>
public Airy1830()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
}

28
GeoUK/Ellipsoids/Airy1830Modified.cs
View File

@@ -1,19 +1,19 @@
namespace GeoUK.Ellipsoids
{
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an Airy1930 ellipsoid and is provided for convienience.
/// </summary>
public class Airy1830Modified : Ellipsoid
{
private const double C_SEMI_MAJOR_AXIS = 6377340.189; //a
private const double C_SEMI_MINOR_AXIS = 6356034.447; //b
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an Airy1930 ellipsoid and is provided for convenience.
/// </summary>
public class Airy1830Modified : Ellipsoid
{
private const double C_SEMI_MAJOR_AXIS = 6377340.189;
private const double C_SEMI_MINOR_AXIS = 6356034.447;
/// <summary>
/// Constructor.
/// </summary>
/// <summary>
/// Constructor.
/// </summary>
public Airy1830Modified()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
}

70
GeoUK/Ellipsoids/Ellipsoid.cs
View File

@@ -7,10 +7,6 @@ namespace GeoUK.Ellipsoids
/// </summary>
public class Ellipsoid
{
private double _semiMajorAxis = 0.0;
private double _semiMinorAxis = 0.0;
private double _eccentricity = 0.0;
private double _eccentricitySquared = 0.0;
/// <summary>
/// Constructor.
@@ -19,10 +15,10 @@ namespace GeoUK.Ellipsoids
/// <param name="semiMinorAxis"></param>
public Ellipsoid(double semiMajorAxis, double semiMinorAxis)
{
_semiMinorAxis = semiMinorAxis;
_semiMajorAxis = semiMajorAxis;
_eccentricitySquared = (Math.Pow(semiMajorAxis, 2) - Math.Pow(semiMinorAxis, 2)) / Math.Pow(semiMajorAxis, 2);
_eccentricity = Math.Sqrt(_eccentricitySquared);
SemiMinorAxis = semiMinorAxis;
SemiMajorAxis = semiMajorAxis;
EccentricitySquared = (Math.Pow(semiMajorAxis, 2) - Math.Pow(semiMinorAxis, 2)) / Math.Pow(semiMajorAxis, 2);
Eccentricity = Math.Sqrt(EccentricitySquared);
}
/// <summary>
@@ -33,88 +29,46 @@ namespace GeoUK.Ellipsoids
public double GetRadiusOfCurvatureInPV(double degreesLatitude)
{
double dblRadians = DegreesToRadians(degreesLatitude);
//this is the working VB example
//GetRadiusOfCurvatureInPV = C_SEMI_MAJOR_AXIS / (1 - C_ECCENTRICITY_SQUARED * Sin(Latitude) ^ 2) ^ 0.5
//GetRadiusOfCurvatureInPV = 6392142.007676
//return (C_SEMI_MAJOR_AXIS / Math.Pow((1 - m_EccentricitySquared * Math.Pow(Math.Sin(dblRadians),2)), 0.5));
return _semiMajorAxis / Math.Pow((1 - _eccentricitySquared * Math.Pow(Math.Sin(dblRadians), 2)), 0.5);
return SemiMajorAxis / Math.Pow((1 - EccentricitySquared * Math.Pow(Math.Sin(dblRadians), 2)), 0.5);
}
/// <summary>
/// Returns the semi-major axis of the ellipsoid.
/// </summary>
public double SemiMajorAxis
{
get
{
return _semiMajorAxis;
}
}
public double SemiMajorAxis { get; }
/// <summary>
/// Returns the semi-major axis of the ellipsoid.
/// </summary>
public double SemiMinorAxis
{
get
{
return _semiMinorAxis;
}
}
public double SemiMinorAxis { get; }
/// <summary>
/// returns the eccentricity of the ellipsoid.
/// </summary>
public double Eccentricity
{
get
{
return _eccentricity;
}
}
public double Eccentricity { get; }
/// <summary>
/// returns the eccentricity squared of the ellipsoid.
/// </summary>
public double EccentricitySquared
{
get
{
return _eccentricitySquared;
}
}
public double EccentricitySquared { get; }
/// <summary>
/// returns the second eccentricity squared of the ellipsoid.
/// </summary>
public double SecondEccentricitySquared
{
get
{
return (Math.Pow(_semiMajorAxis, 2) - Math.Pow(_semiMinorAxis, 2)) / Math.Pow(_semiMinorAxis, 2);
}
}
public double SecondEccentricitySquared => (Math.Pow(SemiMajorAxis, 2) - Math.Pow(SemiMinorAxis, 2)) / Math.Pow(SemiMinorAxis, 2);
/// <summary>
/// Returns radians for a given value of degrees.
/// </summary>
/// <param name="degrees"></param>
/// <returns></returns>
protected double DegreesToRadians(double degrees)
{
return degrees * (Math.PI / 180);
}
protected double DegreesToRadians(double degrees) => degrees * (Math.PI / 180);
/// <summary>
/// Returns degrees for a given value of radians.
/// </summary>
/// <param name="radians"></param>
/// <returns></returns>
protected double RadiansToDegrees(double radians)
{
return radians * (180 / Math.PI);
}
protected double RadiansToDegrees(double radians) => radians * (180 / Math.PI);
}
}

31
GeoUK/Ellipsoids/Grs80.cs
View File

@@ -1,21 +1,16 @@
namespace GeoUK.Ellipsoids
{
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an GRS80 ellipsoid and is provided for convienience.
/// </summary>
public class Grs80 : Ellipsoid
{
//WGS constants
private const double C_SEMI_MAJOR_AXIS = 6378137; //a
private const double C_SEMI_MINOR_AXIS = 6356752.314; //b
/// <summary>
/// Consructor.
/// </summary>
public Grs80()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an GRS80 ellipsoid and is provided for convenience.
/// </summary>
public class Grs80 : Ellipsoid
{
//WGS constants
private const double C_SEMI_MAJOR_AXIS = 6378137;
private const double C_SEMI_MINOR_AXIS = 6356752.314;
public Grs80()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
}

9
GeoUK/Ellipsoids/Hayford1909.cs
View File

@@ -3,13 +3,8 @@
public class Hayford1909 : Ellipsoid
{
//WGS constants
private const double C_SEMI_MAJOR_AXIS = 6378388; //a
private const double C_SEMI_MINOR_AXIS = 6356911.946; //b
/// <summary>
/// Constructor.
/// </summary>
private const double C_SEMI_MAJOR_AXIS = 6378388;
private const double C_SEMI_MINOR_AXIS = 6356911.946;
public Hayford1909()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{

31
GeoUK/Ellipsoids/Wgs84.cs
View File

@@ -1,21 +1,16 @@
namespace GeoUK.Ellipsoids
{
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an WGS84 ellipsoid and is provided for convienience.
/// </summary>
public class Wgs84 : Ellipsoid
{
//WGS constants
private const double C_SEMI_MAJOR_AXIS = 6378137; //a
private const double C_SEMI_MINOR_AXIS = 6356752.3141; //b
/// <summary>
/// Constructor.
/// </summary>
public Wgs84()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
/// <summary>
/// This immutable class, derived from Ellipsoid, represents an WGS84 ellipsoid and is provided for convenience.
/// </summary>
public class Wgs84 : Ellipsoid
{
//WGS constants
private const double C_SEMI_MAJOR_AXIS = 6378137;
private const double C_SEMI_MINOR_AXIS = 6356752.3141;
public Wgs84()
: base(C_SEMI_MAJOR_AXIS, C_SEMI_MINOR_AXIS)
{
}
}
}

5
GeoUK/MathEx.cs
View File

@@ -4,10 +4,7 @@ namespace GeoUK
{
public static class MathEx
{
public static double Secant(double x)
{
return 1 / Math.Cos(x);
}
public static double Secant(double x) => 1 / Math.Cos(x);
// Cosecant Cosec(X) = 1 / Sin(X)
// Cotangent Cotan(X) = 1 / Tan(X)

24
GeoUK/OSTN02Transform.cs
View File

@@ -17,7 +17,7 @@ namespace GeoUK
public static Osgb36 Etrs89ToOsgb(LatitudeLongitude coordinates)
{
EastingNorthing enCoordinates = Convert.ToEastingNorthing(new Grs80(), new BritishNationalGrid(), coordinates);
return Etrs89ToOsgb(enCoordinates, coordinates.ElipsoidalHeight);
return Etrs89ToOsgb(enCoordinates, coordinates.EllipsoidalHeight);
}
private static Osgb36 Etrs89ToOsgb(EastingNorthing coordinates, double ellipsoidHeight)
@@ -48,12 +48,12 @@ namespace GeoUK
string csvRecord = tr.ReadLine();
for (int index = 0; index < 4; index++)
{
if (csvRecord.StartsWith(recordNumbers[index].ToString().Trim() + ",", StringComparison.Ordinal))
{
//dont use add as we need to keep these in same order as record numbers
records[index] = csvRecord;
recordsFound++;
}
if (!csvRecord.StartsWith(recordNumbers[index].ToString().Trim() + ",", StringComparison.Ordinal))
continue;
//don't use add as we need to keep these in same order as record numbers
records[index] = csvRecord;
recordsFound++;
}
}
@@ -103,14 +103,8 @@ namespace GeoUK
/// <param name="eastIndex"></param>
/// <param name="northIndex"></param>
/// <returns></returns>
private static int CalculateRecordNumber(int eastIndex, int northIndex)
{
return eastIndex + (northIndex * 701) + 1;
}
private static int CalculateRecordNumber(int eastIndex, int northIndex) => eastIndex + (northIndex * 701) + 1;
private static Stream GetEmbeddedOSTN02()
{
return ResourceManager.GetEmbeddedResourceStream(typeof(OSTN02Transform).GetTypeInfo().Assembly, "OSTN02_OSGM02_GB.txt");
}
private static Stream GetEmbeddedOSTN02() => ResourceManager.GetEmbeddedResourceStream(typeof(OSTN02Transform).GetTypeInfo().Assembly, "OSTN02_OSGM02_GB.txt");
}
}

5
GeoUK/Projections/BritishNationalGrid.cs
View File

@@ -1,7 +1,7 @@
namespace GeoUK.Projections
{
/// <summary>
/// This immutable class, derived from Projection, represents the British National Grid projection and is provided for convienience.
/// This immutable class, derived from Projection, represents the British National Grid projection and is provided for convenience.
/// </summary>
public class BritishNationalGrid : Projection
{
@@ -11,9 +11,6 @@ namespace GeoUK.Projections
private const double TO_LONG = -2;
private const double N_NH = -100000;
/// <summary>
/// Constructor.
/// </summary>
public BritishNationalGrid()
: base(SCALE_FACTOR, E, N_NH, TO_LAT, TO_LONG)
{

31
GeoUK/Projections/IrishNationalGrid.cs
View File

@@ -1,22 +1,19 @@
namespace GeoUK.Projections
{
/// <summary>
/// This immutable class, derived from Projection, represents the British National Grid projection and is provided for convienience.
/// </summary>
public class IrishNationalGrid : Projection
{
private const double SCALE_FACTOR = 1.000035;
private const double E = 200000;
private const double TO_LAT = 53.5;
private const double TO_LONG = -8;
private const double N_NH = 250000;
/// <summary>
/// This immutable class, derived from Projection, represents the British National Grid projection and is provided for convenience.
/// </summary>
public class IrishNationalGrid : Projection
{
private const double SCALE_FACTOR = 1.000035;
private const double E = 200000;
private const double TO_LAT = 53.5;
private const double TO_LONG = -8;
private const double N_NH = 250000;
/// <summary>
/// Constructor.
/// </summary>
public IrishNationalGrid()
: base(SCALE_FACTOR, E, N_NH, TO_LAT, TO_LONG)
{
}
}
: base(SCALE_FACTOR, E, N_NH, TO_LAT, TO_LONG)
{
}
}
}

73
GeoUK/Projections/Projection.cs
View File

@@ -8,72 +8,40 @@ namespace GeoUK.Projections
/// </summary>
public class Projection
{
private double _trueOriginLatitude = 0.0;
private double _trueOriginLongitude = 0.0;
private double _scaleFactor = 0.0;
private double _trueOriginEasting = 0.0;
private double _trueOriginNorthing = 0.0;
/// <summary>
/// Constructor.
/// </summary>
public Projection(double scaleFactor, double trueOriginEasting, double trueOriginNorthing, double trueOriginLatitude, double trueOriginLongitude)
{
_scaleFactor = scaleFactor;
_trueOriginEasting = trueOriginEasting;
_trueOriginNorthing = trueOriginNorthing;
_trueOriginLatitude = trueOriginLatitude;
_trueOriginLongitude = trueOriginLongitude;
ScaleFactor = scaleFactor;
TrueOriginEasting = trueOriginEasting;
TrueOriginNorthing = trueOriginNorthing;
TrueOriginLatitude = trueOriginLatitude;
TrueOriginLongitude = trueOriginLongitude;
}
/// <summary>
/// Returns the scale factor.
/// </summary>
public double ScaleFactor
{
get { return _scaleFactor; }
}
public double ScaleFactor { get; }
/// <summary>
/// Returns the Easting coordinate of the true origin.
/// </summary>
public double TrueOriginEasting
{
get { return _trueOriginEasting; }
}
public double TrueOriginEasting { get; }
/// <summary>
/// Returns the Northing coordinate of the true origin.
/// </summary>
public double TrueOriginNorthing
{
get { return _trueOriginNorthing; }
}
public double TrueOriginNorthing { get; }
/// <summary>
/// Returns the Latitude coordinate of the true origin for the southern hemisphere.
/// </summary>
public double TrueOriginLatitude
{
get { return _trueOriginLatitude; }
}
public double TrueOriginLatitude { get; }
/// <summary>
/// Returns the Longitude coordinate of the true origin for the southern hemisphere.
/// </summary>
public double TrueOriginLongitude
{
get { return _trueOriginLongitude; }
}
public double TrueOriginLongitude { get; }
//public static double DegreesToRadians(double degrees)
//{
// return degrees * (Math.PI / 180);
//}
//public static double RadiansToDegrees(double radians)
//{
// return radians * (180 / Math.PI);
//}
/// <summary>
/// Returns the integer portion of a division operation.
/// </summary>
@@ -82,32 +50,25 @@ namespace GeoUK.Projections
/// <returns></returns>
protected static double Div(double value, double divisor)
{
double dblResult = 0;
bool blnNegative = false;
//make the division
dblResult = value / divisor;
double dblResult = value / divisor;
//do all calculations on positive numbers
bool blnNegative = false;
if (dblResult < 0)
{
blnNegative = true;
dblResult = dblResult * -1;
dblResult *= -1;
}
//see if there is any remainder
if (dblResult % 1 > 0)
{
dblResult = Math.Ceiling(dblResult) - 1;
}
else
{
dblResult = System.Convert.ToInt32(dblResult, CultureInfo.InvariantCulture);
}
dblResult = dblResult % 1 > 0
? Math.Ceiling(dblResult) - 1
: System.Convert.ToInt32(dblResult, CultureInfo.InvariantCulture);
if (blnNegative)
{
dblResult = dblResult * -1.0;
dblResult *= -1.0;
}
return dblResult;

109
GeoUK/Projections/UniversalTransverseMercator.cs
View File

@@ -2,72 +2,53 @@ using System.Globalization;
namespace GeoUK.Projections
{
/// <summary>
/// This immutable class, derived from Projection, represents the Universal Transverse Mercator projection.
/// The class handles all medidian calculations internally.
/// </summary>
public class UniversalTransverseMercator : Projection
{
private const double SCALE_FACTOR = 0.9996;
private const double E = 500000;
private const double TO_LAT = 0; //is this for britain? These will need to be calculated based on the lat and long being transfformed.
private const double N_NH = 0;
private const double N_SH = 10000000;
/// <summary>
/// This immutable class, derived from Projection, represents the Universal Transverse Mercator projection.
/// The class handles all medidian calculations internally.
/// </summary>
public class UniversalTransverseMercator : Projection
{
private const double SCALE_FACTOR = 0.9996;
private const double E = 500000;
private const double TO_LAT = 0; //is this for Britain? These will need to be calculated based on the lat and long being transformed.
private const double N_NH = 0;
private const double N_SH = 10000000;
/// <summary>
/// Constructor
/// </summary>
/// <param name="degreesLatitude"></param>
/// <param name="degreesLongitude"></param>
public UniversalTransverseMercator(double degreesLatitude, double degreesLongitude)
: base(SCALE_FACTOR, E, CalculateOriginNorthing(degreesLatitude, degreesLongitude), TO_LAT, CalculateLongitudeOrigin(degreesLatitude, degreesLongitude))
{
}
: base(SCALE_FACTOR, E, CalculateOriginNorthing(degreesLatitude), TO_LAT, CalculateLongitudeOrigin(degreesLongitude))
{
}
private static double CalculateOriginNorthing(double degreesLatitude, double degreesLongitude)
{
double falseNorthing = 0.0;
//sort out false northing
if (degreesLatitude < 0)
{
falseNorthing = N_SH;
}
else
{
falseNorthing = N_NH;
}
return falseNorthing;
}
private static double CalculateOriginNorthing(double degreesLatitude) => degreesLatitude < 0 ? N_SH : N_NH;
private static double CalculateLongitudeOrigin(double degreesLatitude, double degreesLongitude)
{
double longOrigin = 0.0;
if (degreesLongitude < 0)
{
if (degreesLongitude > -6)
{
longOrigin = -3;
}
else
{
longOrigin = Div(degreesLongitude, 6.0);
longOrigin = (System.Convert.ToInt32(longOrigin, CultureInfo.InvariantCulture)) * 6 - 3;
}
}
else
{
if (degreesLongitude < 6)
{
longOrigin = 3.0;
}
else
{
longOrigin = Div(degreesLongitude, 6);
longOrigin = (System.Convert.ToInt32(longOrigin, CultureInfo.InvariantCulture)) * 6 + 3;
}
}
return longOrigin;
}
}
private static double CalculateLongitudeOrigin(double degreesLongitude)
{
double longOrigin;
if (degreesLongitude < 0)
{
if (degreesLongitude > -6)
{
longOrigin = -3;
}
else
{
longOrigin = Div(degreesLongitude, 6.0);
longOrigin = (System.Convert.ToInt32(longOrigin, CultureInfo.InvariantCulture)) * 6 - 3;
}
}
else
{
if (degreesLongitude < 6)
{
longOrigin = 3.0;
}
else
{
longOrigin = Div(degreesLongitude, 6);
longOrigin = (System.Convert.ToInt32(longOrigin, CultureInfo.InvariantCulture)) * 6 + 3;
}
}
return longOrigin;
}
}
}

6
GeoUK/ResourceManager.cs
View File

@@ -23,12 +23,12 @@ namespace GeoUK
if (!resourcePaths.Any())
{
throw new Exception(string.Format("Resource ending with {0} not found.", resourceFileName));
throw new Exception($"Resource ending with {resourceFileName} not found.");
}
if (resourcePaths.Count() > 1)
if (resourcePaths.Length > 1)
{
throw new Exception(string.Format("Multiple resources ending with {0} found: {1}{2}", resourceFileName, Environment.NewLine, string.Join(Environment.NewLine, resourcePaths)));
throw new Exception($"Multiple resources ending with {resourceFileName} found: {Environment.NewLine}{string.Join(Environment.NewLine, resourcePaths)}");
}
return assembly.GetManifestResourceStream(resourcePaths.Single());

438
GeoUK/Transform.cs
View File

@@ -3,275 +3,263 @@ using System;
namespace GeoUK
{
/// <summary>
/// This class performs transformations between datums.
/// </summary>
public static class Transform
{
/// <summary>
/// Performs an ETRS89 to OSGB36 datum transformation. Accuracy is approximately 5 meters in all directions.
/// For this method, ERTS89, ITRS2000 and WGS84 datums can be considered the same.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <remarks>
/// This method uses a Helmert transformation to determine the OSGB36 coordinates.
/// Whilst only accurate to 5 meters in all directions, it is extremely fast.
/// </remarks>
/// <summary>
/// This class performs transformations between datums.
/// </summary>
public static class Transform
{
/// <summary>
/// Performs an ETRS89 to OSGB36 datum transformation. Accuracy is approximately 5 meters in all directions.
/// For this method, ERTS89, ITRS2000 and WGS84 datums can be considered the same.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <remarks>
/// This method uses a Helmert transformation to determine the OSGB36 coordinates.
/// Whilst only accurate to 5 meters in all directions, it is extremely fast.
/// </remarks>
public static Cartesian Etrs89ToOsgb36(Cartesian coordinates)
{
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//-446.448 125.157 - 542.060 20.4894 - 0.1502 - 0.247 - 0.8421
{
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//-446.448 125.157 - 542.060 20.4894 - 0.1502 - 0.247 - 0.8421
//set up the parameters
double tx = -446.448;
double ty = 125.157;
double tz = -542.060;
double s = 20.4894;
double rx = ToRadians(ToDecimelDegrees(0, 0, -0.1502));
double ry = ToRadians(ToDecimelDegrees(0, 0, -0.247));
double rz = ToRadians(ToDecimelDegrees(0, 0, -0.8421));
//set up the parameters
double tx = -446.448;
double ty = 125.157;
double tz = -542.060;
double s = 20.4894;
double rx = ToRadians(ToDecimalDegrees(0, 0, -0.1502));
double ry = ToRadians(ToDecimalDegrees(0, 0, -0.247));
double rz = ToRadians(ToDecimalDegrees(0, 0, -0.8421));
Cartesian result = HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
Cartesian result = HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
return result;
}
return result;
}
/// <summary>
/// <summary>
/// Performs an OSGB36 to ETRS89 datum transformation. Accuracy is approximately 5 meters in all directions.
/// For this method, ERTS89, ITRS2000 and WGS84 datums can be considered the same.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <returns></returns>
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <returns></returns>
public static Cartesian Osgb36ToEtrs89(Cartesian coordinates)
{
//(BUT CHANGE SIGNS OF EACH PARAMETER FOR REVERSE)
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//-446.448 125.157 - 542.060 20.4894 - 0.1502 - 0.247 - 0.8421
{
//(BUT CHANGE SIGNS OF EACH PARAMETER FOR REVERSE)
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//-446.448 125.157 - 542.060 20.4894 - 0.1502 - 0.247 - 0.8421
double tx = 446.448;
double ty = -125.157;
double tz = 542.060;
double s = -20.4894;
double rx = ToRadians(ToDecimelDegrees(0, 0, 0.1502));
double ry = ToRadians(ToDecimelDegrees(0, 0, 0.247));
double rz = ToRadians(ToDecimelDegrees(0, 0, 0.8421));
double tx = 446.448;
double ty = -125.157;
double tz = 542.060;
double s = -20.4894;
double rx = ToRadians(ToDecimalDegrees(0, 0, 0.1502));
double ry = ToRadians(ToDecimalDegrees(0, 0, 0.247));
double rz = ToRadians(ToDecimalDegrees(0, 0, 0.8421));
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
/*
/*
* //(BUT CHANGE SIGNS OF EACH PARAMETER FOR REVERSE)
tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
-446.448 +125.157 -542.060 +20.4894 -0.1502 -0.2470 -0.8421
*/
/// <summary>
/// Performs an ETRS89 to ITRS2000 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Etrs89ToItrs2000(Cartesian coordinates, int epochYear)
{
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.054 0.051 - 0.048 0 0.000081 dt 0.00049 dt - 0.000792 dt
/// <summary>
/// Performs an ETRS89 to ITRS2000 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Etrs89ToItrs2000(Cartesian coordinates, int epochYear)
{
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.054 0.051 - 0.048 0 0.000081 dt 0.00049 dt - 0.000792 dt
//dt represents shift in years since time of survey to when ETRS89 was determined
//dt represents shift in years since time of survey to when ETRS89 was determined
int dt = epochYear - 1989;
//set up the parameters
double tx = Negate(0.054);
double ty = Negate(0.051);
double tz = Negate(-0.048);
double s = 0;
double rx = Negate(ToRadians(ToDecimelDegrees(0, 0, 0.000081) * dt));
double ry = Negate(ToRadians(ToDecimelDegrees(0, 0, 0.00049) * dt));
double rz = Negate(ToRadians(ToDecimelDegrees(0, 0, -0.000792) * dt));
int dt = epochYear - 1989;
//set up the parameters
double tx = Negate(0.054);
double ty = Negate(0.051);
double tz = Negate(-0.048);
double s = 0;
double rx = Negate(ToRadians(ToDecimalDegrees(0, 0, 0.000081) * dt));
double ry = Negate(ToRadians(ToDecimalDegrees(0, 0, 0.00049) * dt));
double rz = Negate(ToRadians(ToDecimalDegrees(0, 0, -0.000792) * dt));
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
/// <summary>
/// Performs an ITRS2000 to ETRS89 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Itrs2000ToEtrs89(Cartesian coordinates, int epochYear)
{
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.054 0.051 - 0.048 0 0.000081 dt 0.00049 dt - 0.000792 dt
/// <summary>
/// Performs an ITRS2000 to ETRS89 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Itrs2000ToEtrs89(Cartesian coordinates, int epochYear)
{
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.054 0.051 - 0.048 0 0.000081 dt 0.00049 dt - 0.000792 dt
//dt represents shift in years since time of survey to when ETRS89 was determined
//dt represents shift in years since time of survey to when ETRS89 was determined
int dt = epochYear - 1989;
//set up the parameters
double tx = 0.054;
double ty = 0.051;
double tz = -0.048;
double s = 0;
double rx = ToRadians(ToDecimelDegrees(0, 0, 0.000081) * dt);
double ry = ToRadians(ToDecimelDegrees(0, 0, 0.00049) * dt);
double rz = ToRadians(ToDecimelDegrees(0, 0, -0.000792) * dt);
int dt = epochYear - 1989;
//set up the parameters
double tx = 0.054;
double ty = 0.051;
double tz = -0.048;
double s = 0;
double rx = ToRadians(ToDecimalDegrees(0, 0, 0.000081) * dt);
double ry = ToRadians(ToDecimalDegrees(0, 0, 0.00049) * dt);
double rz = ToRadians(ToDecimalDegrees(0, 0, -0.000792) * dt);
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
/// <summary>
/// Performs an ITRS94/96/97 to ETRS89 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Itrs97ToEtrs89(Cartesian coordinates, int epochYear)
{
//ITRS94/96/97 to ETRS89 datum transformation
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.041 0.041 - 0.049 0 0.00020 dt 0.00050 dt - 0.00065 dt
/// <summary>
/// Performs an ITRS94/96/97 to ETRS89 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Itrs97ToEtrs89(Cartesian coordinates, int epochYear)
{
//ITRS94/96/97 to ETRS89 datum transformation
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.041 0.041 - 0.049 0 0.00020 dt 0.00050 dt - 0.00065 dt
//dt represents shift in years since time of survey to when ETRS89 was determined
//dt represents shift in years since time of survey to when ETRS89 was determined
int dt = epochYear - 1989;
//set up the parameters
double tx = 0.041;
double ty = 0.041;
double tz = -0.049;
double s = 0;
double rx = ToRadians(ToDecimelDegrees(0, 0, 0.00020) * dt);
double ry = ToRadians(ToDecimelDegrees(0, 0, 0.00050) * dt);
double rz = ToRadians(ToDecimelDegrees(0, 0, -0.00065) * dt);
int dt = epochYear - 1989;
//set up the parameters
double tx = 0.041;
double ty = 0.041;
double tz = -0.049;
double s = 0;
double rx = ToRadians(ToDecimalDegrees(0, 0, 0.00020) * dt);
double ry = ToRadians(ToDecimalDegrees(0, 0, 0.00050) * dt);
double rz = ToRadians(ToDecimalDegrees(0, 0, -0.00065) * dt);
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
/// <summary>
/// Performs an ETRS89 to ITRS94/96/97 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Etrs89ToItrs97(Cartesian coordinates, int epochYear)
{
//ITRS94/96/97 to ETRS89 datum transformation (BUT CHANGE SIGNS OF EACH PARAMETER FOR REVERSE)
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.041 0.041 - 0.049 0 0.00020 dt 0.00050 dt - 0.00065 dt
/// <summary>
/// Performs an ETRS89 to ITRS94/96/97 datum transformation.
/// </summary>
/// <param name="coordinates">Cartesian Coordinates to be transformed.</param>
/// <param name="epochYear">Refers to the year the data specified in coordinates was gathered.</param>
/// <returns></returns>
public static Cartesian Etrs89ToItrs97(Cartesian coordinates, int epochYear)
{
//ITRS94/96/97 to ETRS89 datum transformation (BUT CHANGE SIGNS OF EACH PARAMETER FOR REVERSE)
//tX (m) tY (m) tZ (m) s (ppm) rX (sec) rY (sec) rZ (sec)
//0.041 0.041 - 0.049 0 0.00020 dt 0.00050 dt - 0.00065 dt
//dt represents shift in years since time of survey to when ETRS89 was determined
//dt represents shift in years since time of survey to when ETRS89 was determined
int dt = epochYear - 1989;
//set up the parameters
double tx = Negate(0.041);
double ty = Negate(0.041);
double tz = Negate(-0.049);
double s = 0;
double rx = Negate(ToRadians(ToDecimelDegrees(0, 0, 0.00020) * dt));
double ry = Negate(ToRadians(ToDecimelDegrees(0, 0, 0.00050) * dt));
double rz = Negate(ToRadians(ToDecimelDegrees(0, 0, -0.00065) * dt));
int dt = epochYear - 1989;
//set up the parameters
double tx = Negate(0.041);
double ty = Negate(0.041);
double tz = Negate(-0.049);
double s = 0;
double rx = Negate(ToRadians(ToDecimalDegrees(0, 0, 0.00020) * dt));
double ry = Negate(ToRadians(ToDecimalDegrees(0, 0, 0.00050) * dt));
double rz = Negate(ToRadians(ToDecimalDegrees(0, 0, -0.00065) * dt));
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
return HelmertTransformation(coordinates, tx, ty, tz, rx, ry, rz, s);
}
private static double ToDecimelDegrees(int degrees, int minutes, double seconds)
{
//determine seconds as minutes
double m = minutes + (seconds / 60.0);
return ToDecimelDegrees(degrees, m);
}
private static double ToDecimalDegrees(int degrees, int minutes, double seconds)
{
//determine seconds as minutes
double m = minutes + (seconds / 60.0);
return ToDecimalDegrees(degrees, m);
}
private static double ToDecimelDegrees(int degrees, double minutes)
{
//determine minutes as derees
return degrees + (minutes / 60.0);
}
// determine minutes as degrees
private static double ToDecimalDegrees(int degrees, double minutes) => degrees + (minutes / 60.0);
private static double ToRadians(double degrees)
{
return degrees * (Math.PI / 180.0);
}
private static double ToRadians(double degrees) => degrees * (Math.PI / 180.0);
private static double ToDegrees(double radians)
{
return radians * (180.0 / Math.PI);
}
private static double ToDegrees(double radians) => radians * (180.0 / Math.PI);
/// <summary>
/// This seven parameter method can be used to transform coordinates between datums.
/// </summary>
/// <remarks>
/// This method assumes that the rotation
/// parameters are <20>small<6C>. Rotation parameters between geodetic cartesian systems are usually less than 5
/// seconds of arc, because the axes are conventionally aligned to the Greenwich Meridian and the Pole.
/// Do not use this formula for larger angles.
/// </remarks>
/// <param name="coordinates"></param>
/// <param name="translationX"></param>
/// <param name="translationY"></param>
/// <param name="translationZ"></param>
/// <param name="rotationX"></param>
/// <param name="rotationY"></param>
/// <param name="rotationZ"></param>
/// <param name="scaleFactorPpm"></param>
/// <returns></returns>
public static Cartesian HelmertTransformation(Cartesian coordinates, double translationX, double translationY, double translationZ, double rotationX, double rotationY, double rotationZ, double scaleFactorPpm)
{
//to compute this helmert translation we have to multiply XYZa by R and add to T
//
// : X :B : TX : : 1+s -rz ry : : X :A
// : : : : : : : :
// : Y : = : TY : + : rz 1+s -rx :. : Y :
// : : : : : : : :
// : Z : : TZ : : -ry rx 1+s : : Z :
/// <summary>
/// This seven parameter method can be used to transform coordinates between datums.
/// </summary>
/// <remarks>
/// This method assumes that the rotation
/// parameters are <20>small<6C>. Rotation parameters between geodetic cartesian systems are usually less than 5
/// seconds of arc, because the axes are conventionally aligned to the Greenwich Meridian and the Pole.
/// Do not use this formula for larger angles.
/// </remarks>
/// <param name="coordinates"></param>
/// <param name="translationX"></param>
/// <param name="translationY"></param>
/// <param name="translationZ"></param>
/// <param name="rotationX"></param>
/// <param name="rotationY"></param>
/// <param name="rotationZ"></param>
/// <param name="scaleFactorPpm"></param>
/// <returns></returns>
public static Cartesian HelmertTransformation(Cartesian coordinates, double translationX, double translationY, double translationZ, double rotationX, double rotationY, double rotationZ, double scaleFactorPpm)
{
//to compute this helmert translation we have to multiply XYZa by R and add to T
//
// : X :B : TX : : 1+s -rz ry : : X :A
// : : : : : : : :
// : Y : = : TY : + : rz 1+s -rx :. : Y :
// : : : : : : : :
// : Z : : TZ : : -ry rx 1+s : : Z :
//scale factor passed in as a parts/million measure
double scaleFactor = scaleFactorPpm / 1000000.0;
//scale factor passed in as a parts/million measure
double scaleFactor = scaleFactorPpm / 1000000.0;
//create initial matrixes (cols, rows)
double[] XYZa = new double[3]; //initial xyz
XYZa[0] = coordinates.X;
XYZa[1] = coordinates.Y;
XYZa[2] = coordinates.Z;
//create initial matrixes (cols, rows)
double[] XYZa = new double[3]; //initial xyz
XYZa[0] = coordinates.X;
XYZa[1] = coordinates.Y;
XYZa[2] = coordinates.Z;
//populate main matrix 'R'
double[,] R = new double[3, 3];
//top row
R[0, 0] = 1 + scaleFactor; R[1, 0] = -rotationZ; R[2, 0] = rotationY;
//second row
R[0, 1] = rotationZ; R[1, 1] = 1 + scaleFactor; R[2, 1] = -rotationX;
//third row
R[0, 2] = -rotationY; R[1, 2] = rotationX; R[2, 2] = 1 + scaleFactor;
//populate main matrix 'R'
double[,] R = new double[3, 3];
//top row
R[0, 0] = 1 + scaleFactor; R[1, 0] = -rotationZ; R[2, 0] = rotationY;
//second row
R[0, 1] = rotationZ; R[1, 1] = 1 + scaleFactor; R[2, 1] = -rotationX;
//third row
R[0, 2] = -rotationY; R[1, 2] = rotationX; R[2, 2] = 1 + scaleFactor;
//populate matrix 'T'
double[] T = new double[3];
T[0] = translationX; T[1] = translationY; T[2] = translationZ;
//populate matrix 'T'
double[] T = new double[3];
T[0] = translationX; T[1] = translationY; T[2] = translationZ;
//intermediate result of the multiplication goes here
double[] temp = new double[3];
//intermediate result of the multiplication goes here
double[] temp = new double[3];
//final result goes here
//double[,] XYZb = new double[1, 3]; //result
//final result goes here
//double[,] XYZb = new double[1, 3]; //result
//start with the multiplication of matrix XYZa and R
for (int row = 0; row < 3; row++)
{
for (int col = 0; col < 3; col++)
{
temp[row] = temp[row] + (R[col, row] * XYZa[col]);
}
}
//start with the multiplication of matrix XYZa and R
for (int row = 0; row < 3; row++)
{
for (int col = 0; col < 3; col++)
{
temp[row] = temp[row] + (R[col, row] * XYZa[col]);
}
}
//adding to T whilst creating the cartesian coordinates
return new Cartesian(T[0] + temp[0], T[1] + temp[1], T[2] + temp[2]);
}
//adding to T whilst creating the cartesian coordinates
return new Cartesian(T[0] + temp[0], T[1] + temp[1], T[2] + temp[2]);
}
/// <summary>
/// Helper funtion to reverse the sign of a value. Helps code to be more readable.
/// </summary>
/// <param name="value"></param>
/// <returns></returns>
private static double Negate(double value)
{
return value * -1.0;
}
}
/// <summary>
/// Helper function to reverse the sign of a value. Helps code to be more readable.
/// </summary>
/// <param name="value"></param>
/// <returns></returns>
private static double Negate(double value) => value * -1.0;
}
}