Updated azimuthal projections

This commit is contained in:
ClemensFischer 2026-01-21 12:08:12 +01:00
parent d5ee7bbfc0
commit 1bf4cab071
6 changed files with 93 additions and 39 deletions

View file

@ -29,14 +29,15 @@ namespace MapControl
public override Matrix RelativeScale(double latitude, double longitude)
{
var p = GetProjectedPoint(latitude, longitude);
var k = 1d;
if (p.CosC < 1d)
if (Math.Abs(p.CosC) >= 1d)
{
var c = Math.Acos(p.CosC);
k = c / Math.Sin(c); // p.195 (25-2)
return new Matrix(1d, 0d, 0d, 1d, 0d, 0d);
}
var c = Math.Acos(p.CosC);
var k = c / Math.Sin(c); // p.195 (25-2)
var scale = new Matrix(1d, 0d, 0d, k, 0d, 0d); // h == 1
scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
return scale;
@ -45,14 +46,15 @@ namespace MapControl
public override Point? LocationToMap(double latitude, double longitude)
{
var p = GetProjectedPoint(latitude, longitude);
var k = 1d;
if (p.CosC < 1d)
if (Math.Abs(p.CosC) >= 1d)
{
var c = Math.Acos(p.CosC);
k = c / Math.Sin(c); // p.195 (25-2)
return new Point();
}
var c = Math.Acos(p.CosC);
var k = c / Math.Sin(c); // p.195 (25-2)
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.195 (22-4/5)
}

View file

@ -46,14 +46,39 @@ namespace MapControl
protected Location GetLocation(double x, double y, double rho, double sinC)
{
var cosC = Math.Sqrt(1d - sinC * sinC);
var cos2C = 1d - sinC * sinC;
if (cos2C < 0d)
{
return null;
}
var cosC = Math.Sqrt(cos2C);
var phi1 = Center.Latitude * Math.PI / 180d;
var cosPhi1 = Math.Cos(phi1);
var sinPhi1 = Math.Sin(phi1);
var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / rho); // (20-14)
var dLambda = Math.Atan2(x * sinC, rho * cosPhi1 * cosC - y * sinPhi1 * sinC); // (20-15), λ - λ0
double u, v;
return new Location(phi * 180d / Math.PI, dLambda * 180d / Math.PI + Center.Longitude);
if (Center.Latitude == 90d) // (20-16)
{
u = x;
v = -y;
}
else if (Center.Latitude == -90d) // (20-17)
{
u = x;
v = y;
}
else // (20-15)
{
u = x * sinC;
v = rho * cosPhi1 * cosC - y * sinPhi1 * sinC;
}
return new Location(
phi * 180d / Math.PI,
Math.Atan2(u, v) * 180d / Math.PI + Center.Longitude);
}
}
}

View file

@ -41,7 +41,7 @@ namespace MapControl
{
var p = GetProjectedPoint(latitude, longitude);
if (p.CosC <= 0d)
if (p.CosC <= 0d) // p.167 "If cos c is zero or negative, the point is to be rejected."
{
return null;
}

View file

@ -6,7 +6,7 @@ namespace MapControl
{
public virtual MapProjection GetProjection(string crsId)
{
MapProjection projection = crsId switch
var projection = crsId switch
{
WebMercatorProjection.DefaultCrsId => new WebMercatorProjection(),
WorldMercatorProjection.DefaultCrsId => new WorldMercatorProjection(),
@ -26,17 +26,27 @@ namespace MapControl
return projection ?? throw new NotSupportedException($"MapProjection \"{crsId}\" is not supported.");
}
public MapProjection GetProjectionFromEpsgCode(string crsId) =>
crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode) ? GetProjection(epsgCode) : null;
public virtual MapProjection GetProjection(int epsgCode) => epsgCode switch
public virtual MapProjection GetProjection(int epsgCode)
{
var code when code >= Etrs89UtmProjection.FirstZoneEpsgCode && code <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(epsgCode % 100),
var code when code >= Nad27UtmProjection.FirstZoneEpsgCode && code <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(epsgCode % 100),
var code when code >= Nad83UtmProjection.FirstZoneEpsgCode && code <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(epsgCode % 100),
var code when code >= Wgs84UtmProjection.FirstZoneNorthEpsgCode && code <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.North),
var code when code >= Wgs84UtmProjection.FirstZoneSouthEpsgCode && code <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.South),
_ => null
};
return epsgCode switch
{
var c when c is >= Etrs89UtmProjection.FirstZoneEpsgCode
and <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(c % 100),
var c when c is >= Nad27UtmProjection.FirstZoneEpsgCode
and <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(c % 100),
var c when c is >= Nad83UtmProjection.FirstZoneEpsgCode
and <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(c % 100),
var c when c is >= Wgs84UtmProjection.FirstZoneNorthEpsgCode
and <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.North),
var c when c is >= Wgs84UtmProjection.FirstZoneSouthEpsgCode
and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South),
_ => null
};
}
protected MapProjection GetProjectionFromEpsgCode(string crsId)
{
return crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode) ? GetProjection(epsgCode) : null;
}
}
}

View file

@ -40,7 +40,12 @@ namespace MapControl
{
var p = GetProjectedPoint(latitude, longitude);
return p.CosC >= 0d ? new Point(EarthRadius * p.X, EarthRadius * p.Y) : null; // p.149 (20-3/4)
if (p.CosC < 0d) // p.149 "If cos c is negative, the point is not to be plotted."
{
return null;
}
return new Point(EarthRadius * p.X, EarthRadius * p.Y); // p.149 (20-3/4)
}
public override Location MapToLocation(double x, double y)
@ -48,7 +53,7 @@ namespace MapControl
var rho = Math.Sqrt(x * x + y * y);
var sinC = rho / EarthRadius; // p.150 (20-19)
return sinC <= 1d ? GetLocation(x, y, rho, sinC) : null;
return GetLocation(x, y, rho, sinC);
}
}
}

View file

@ -27,7 +27,7 @@ namespace MapControl.Projections
public override MapProjection GetProjection(string crsId)
{
MapProjection projection = crsId switch
return crsId switch
{
MapControl.WebMercatorProjection.DefaultCrsId => new WebMercatorProjection(),
MapControl.WorldMercatorProjection.DefaultCrsId => new WorldMercatorProjection(),
@ -36,21 +36,33 @@ namespace MapControl.Projections
MapControl.Wgs84AutoUtmProjection.DefaultCrsId => new Wgs84AutoUtmProjection(),
MapControl.OrthographicProjection.DefaultCrsId => new Wgs84OrthographicProjection(),
MapControl.StereographicProjection.DefaultCrsId => new Wgs84StereographicProjection(),
_ => GetProjectionFromEpsgCode(crsId)
_ => GetProjectionFromEpsgCode(crsId) ?? base.GetProjection(crsId)
};
return projection ?? base.GetProjection(crsId);
}
public override MapProjection GetProjection(int epsgCode) => epsgCode switch
public override MapProjection GetProjection(int epsgCode)
{
var code when code >= Ed50UtmProjection.FirstZoneEpsgCode && code <= Ed50UtmProjection.LastZoneEpsgCode => new Ed50UtmProjection(epsgCode % 100),
var code when code >= Etrs89UtmProjection.FirstZoneEpsgCode && code <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(epsgCode % 100),
var code when code >= Nad27UtmProjection.FirstZoneEpsgCode && code <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(epsgCode % 100),
var code when code >= Nad83UtmProjection.FirstZoneEpsgCode && code <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(epsgCode % 100),
var code when code >= Wgs84UtmProjection.FirstZoneNorthEpsgCode && code <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.North),
var code when code >= Wgs84UtmProjection.FirstZoneSouthEpsgCode && code <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.South),
_ => CoordinateSystemWkts.TryGetValue(epsgCode, out string wkt) ? new ProjNetMapProjection(wkt) : base.GetProjection(epsgCode)
};
if (CoordinateSystemWkts.TryGetValue(epsgCode, out string wkt))
{
return new ProjNetMapProjection(wkt);
}
return epsgCode switch
{
var c when c is >= Ed50UtmProjection.FirstZoneEpsgCode
and <= Ed50UtmProjection.LastZoneEpsgCode => new Ed50UtmProjection(c % 100),
var c when c is >= Etrs89UtmProjection.FirstZoneEpsgCode
and <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(c % 100),
var c when c is >= Nad27UtmProjection.FirstZoneEpsgCode
and <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(c % 100),
var c when c is >= Nad83UtmProjection.FirstZoneEpsgCode
and <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(c % 100),
var c when c is >= Wgs84UtmProjection.FirstZoneNorthEpsgCode
and <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.North),
var c when c is >= Wgs84UtmProjection.FirstZoneSouthEpsgCode
and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South),
_ => base.GetProjection(epsgCode)
};
}
}
}