Skip to content

Instantly share code, notes, and snippets.

@oliverheilig
Last active August 15, 2023 12:41
Show Gist options
  • Save oliverheilig/7029947 to your computer and use it in GitHub Desktop.
Save oliverheilig/7029947 to your computer and use it in GitHub Desktop.
GeoTransform
// This snippet contains methods to transform between the
// various "coordinate formats" used by PTV components
using System;
public class Program
{
public static void Main()
{
var pWgs = new Point(8.4, 49.0); // Karlsruhe
var pPtv_Geodecimal =
GeoTransform.Trans(CoordinateFormat.Wgs84,
CoordinateFormat.Ptv_Geodecimal,
pWgs);
var pPtv_Mercator =
GeoTransform.Trans(CoordinateFormat.Wgs84,
CoordinateFormat.Ptv_Mercator,
pWgs);
var pWeb_Mercator =
GeoTransform.Trans(CoordinateFormat.Wgs84,
CoordinateFormat.Web_Mercator,
pWgs);
var pPtv_SmartUnits =
GeoTransform.Trans(CoordinateFormat.Wgs84,
CoordinateFormat.Ptv_SmartUnits,
pWgs);
Console.WriteLine("Wgs84: {0}/{1}",
pWgs.X, pWgs.Y);
Console.WriteLine("PTV Geodecimal: {0}/{1}",
pPtv_Geodecimal.X, pPtv_Geodecimal.Y);
Console.WriteLine("PTV Mercator: {0}/{1}",
pPtv_Mercator.X, pPtv_Mercator.Y);
Console.WriteLine("Google Mercator: {0}/{1}",
pWeb_Mercator.X, pWeb_Mercator.Y);
Console.WriteLine("PTV Smart Units: {0}/{1}",
pPtv_SmartUnits.X, pPtv_SmartUnits.Y);
}
}
public enum CoordinateFormat
{
/// <summary>
/// Lat/Lon EPGS:4326
/// </summary>
Wgs84,
/// <summary>
/// PTV Spherical Mercator EPSG:505456
/// </summary>
Ptv_Mercator,
/// <summary>
/// PTV Geodecimal (WGS * 100000)
/// </summary>
Ptv_Geodecimal,
/// <summary>
/// PTV Smart Units
/// </summary>
Ptv_SmartUnits,
/// <summary>
/// Spherical Web Mercator Projection
/// (aka Google Mercator) EPSG:900913
/// </summary>
Web_Mercator,
}
public static class GeoTransform
{
public static Point Trans(CoordinateFormat inFormat, CoordinateFormat outFormat, Point point)
{
// int == out -> return
if (inFormat == outFormat)
return point;
// direct transformations
if (inFormat == CoordinateFormat.Ptv_Geodecimal && outFormat == CoordinateFormat.Wgs84)
return Geodecimal_2_WGS84(point);
if (inFormat == CoordinateFormat.Wgs84 && outFormat == CoordinateFormat.Ptv_Geodecimal)
return WGS84_2_Geodecimal(point);
if (inFormat == CoordinateFormat.Ptv_Mercator && outFormat == CoordinateFormat.Ptv_SmartUnits)
return Mercator_2_SmartUnits(point);
if (inFormat == CoordinateFormat.Ptv_SmartUnits && outFormat == CoordinateFormat.Ptv_Mercator)
return SmartUnits_2_Mercator(point);
if (inFormat == CoordinateFormat.Ptv_Mercator && outFormat == CoordinateFormat.Wgs84)
return SphereMercator_2_Wgs(point, Ptv_Radius);
if (inFormat == CoordinateFormat.Wgs84 && outFormat == CoordinateFormat.Ptv_Mercator)
return Wgs_2_SphereMercator(point, Ptv_Radius);
if (inFormat == CoordinateFormat.Web_Mercator && outFormat == CoordinateFormat.Wgs84)
return SphereMercator_2_Wgs(point, Google_Radius);
if (inFormat == CoordinateFormat.Wgs84 && outFormat == CoordinateFormat.Web_Mercator)
return Wgs_2_SphereMercator(point, Google_Radius);
if (inFormat == CoordinateFormat.Ptv_Mercator && outFormat == CoordinateFormat.Web_Mercator)
return Ptv_2_Google(point);
if (inFormat == CoordinateFormat.Web_Mercator && outFormat == CoordinateFormat.Ptv_Mercator)
return Google_2_Ptv(point);
// transitive transformations
if (inFormat == CoordinateFormat.Ptv_SmartUnits)
return Trans(CoordinateFormat.Ptv_Mercator, outFormat, SmartUnits_2_Mercator(point));
if (outFormat == CoordinateFormat.Ptv_SmartUnits)
return Mercator_2_SmartUnits(Trans(inFormat, CoordinateFormat.Ptv_Mercator, point));
if (inFormat == CoordinateFormat.Ptv_Geodecimal)
return Trans(CoordinateFormat.Wgs84, outFormat, Geodecimal_2_WGS84(point));
if (outFormat == CoordinateFormat.Ptv_Geodecimal)
return WGS84_2_Geodecimal(Trans(inFormat, CoordinateFormat.Wgs84, point));
// this should not happen
throw new NotImplementedException(string.Format("transformation not implemented for {0} to {1}",
inFormat.ToString(), outFormat.ToString()));
}
#region geographic formats
public static Point Geodecimal_2_WGS84(Point point)
{
return new Point(point.X / 100000.0, point.Y / 100000.0);
}
public static Point WGS84_2_Geodecimal(Point point)
{
return new Point(point.X * 100000.0, point.Y * 100000.0);
}
#endregion
#region smart units
const double SMARTFACTOR = 4.809543;
const int SMART_OFFSET = 20015087;
public static Point SmartUnits_2_Mercator(Point point)
{
return new Point(
(point.X * SMARTFACTOR) - SMART_OFFSET,
(point.Y * SMARTFACTOR) - SMART_OFFSET);
}
public static Point Mercator_2_SmartUnits(Point point)
{
return new Point(
(point.X + SMART_OFFSET) / SMARTFACTOR,
(point.Y + SMART_OFFSET) / SMARTFACTOR);
}
#endregion
#region Spherical Mercator
const double Ptv_Radius = 6371000.0;
const double Google_Radius = 6378137.0;
public static Point Ptv_2_Google(Point point)
{
return new Point(point.X / Ptv_Radius * Google_Radius, point.Y / Ptv_Radius * Google_Radius);
}
public static Point Google_2_Ptv(Point point)
{
return new Point(point.X / Google_Radius * Ptv_Radius, point.Y / Google_Radius * Ptv_Radius);
}
public static Point Wgs_2_SphereMercator(Point point, double earthRadius)
{
double x = earthRadius * point.X * Math.PI / 180.0;
double y = earthRadius * Math.Log(Math.Tan(Math.PI / 4.0 + point.Y * Math.PI / 360.0));
return new Point(x, y);
}
public static Point SphereMercator_2_Wgs(Point point, double earthRadius)
{
double x = (180.0 / Math.PI) * (point.X / earthRadius);
double y = (360 / Math.PI) * (Math.Atan(Math.Exp(point.Y / earthRadius)) - (Math.PI / 4));
return new Point(x, y);
}
#endregion
}
public class Point
{
public Point(double x, double y)
{
this.X = x;
this.Y = y;
}
public double X { get; set; }
public double Y { get; set; }
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment