How to get the SpatialReference horizon polygon?

[AlanStewart]

There doesn't seem to be a direct way to get this as was possible with ArcObjects. It looks like it should be doable with GeometryEngine.GeodesicEllipse() with the parameters setup correctly. I've done it successfully for one SpatialReference, but so far unsuccessful with others. I need this to be general for all projected coordinate systems, any part of the Earth. The code below worked for NAD 1983 (2011) Contiguous  USA Albers. Any suggestions on how to improve this?

public static Polygon GetSpatialReferenceHorizon(SpatialReference sr)
{
   GeodesicEllipseParameter param = new GeodesicEllipseParameter();

   param.LinearUnit = sr.Unit as LinearUnit;

   param.SemiAxis1Length = 50000000 * sr.Unit.ConversionFactor; // larger than the circumference of the Earth
   param.SemiAxis2Length = 50000000 * sr.Unit.ConversionFactor;
   param.AxisDirection = 0;
   param.Center = new Coordinate2D(0, 0);

   param.VertexCount = 1000;
   param.OutGeometryType = GeometryType.Polygon;

   Polygon horizon = GeometryEngine.Instance.GeodesicEllipse(param, sr) as Polygon;
   return horizon;
}

[AlanStewart]

I still struggle with this. You would think constructing a circle that encompasses the SpatialReference.Domain.Extent would always work, but I generally get an invalid argument exception.

// 'sr' is the passed in SpatialReference

double x = sr.Domain.Extent.Width / 2;
double y = sr.Domain.Extent.Height / 2;
double radius = Math.Sqrt((x * x) + (y * y));

GeodesicEllipseParameter param = new GeodesicEllipseParameter();

param.LinearUnit = sr.Unit as LinearUnit;

param.SemiAxis1Length = radius;
param.SemiAxis2Length = radius;
param.AxisDirection = 0;
param.Center = sr.Domain.Extent.CenterCoordinate;

param.VertexCount = 1000;
param.OutGeometryType = GeometryType.Polygon;

Polygon horizon = GeometryEngine.Instance.GeodesicEllipse(param,
sr)
as Polygon;

[AlanStewart]

*** ESRI BUG ***

For the NAD 1983 (2011) Contiguous  USA Albers spatial reference domain is returned as:

-16901100.0, -6972200.0 and 900703024374.099, 900712953274.099 (min and max points).

However, the horizon polygon returned by GeodesicEllipse(), which appears to be correct, has an extent of -16900932.371123984, -6972041.5589691568 and 16900932.361594178, 15297987.438447423

Could be why I always get an invalid argument exception when I construct an ellipse from the domain.

I'm using Pro 2.4.

[AnnetteLocke]

Hi Alan,

I will take a look and get back to you. 

Thanks,

Annette

[AnnetteLocke]

Alan,

You are right. This won't work for the reason that you stated. We will add a method to get the horizon from the spatial reference in a near-term release. 

Thank you,

Annette

[AlanStewart]

Here is an updated algorithm that does not rely on the domain, but for some reason this only works with 2331 of the spatial references I have tested. The other 613 fail with the std::invalid_argument exception. I do not see what is shared by the ones that fail that is not shared with ones that succeed.

           // "sr" is the passed in SpatialReference

            double a = sr.Datum.SpheroidSemiMajorAxis;
            double b = sr.Datum.SpheroidSemiMinorAxis;
            double circumference = Math.PI * (3 * (+b) - Math.Sqrt((3 * a + b) * (a + 3 * b)));

            GeodesicEllipseParameter param = new GeodesicEllipseParameter();

            param.LinearUnit = sr.Unit as LinearUnit;

            param.SemiAxis1Length = circumference / sr.Unit.ConversionFactor;
            param.SemiAxis2Length = circumference / sr.Unit.ConversionFactor;
            param.AxisDirection = 0;
            param.Center = new Coordinate2D(0, 0);

            param.VertexCount = 1000;
            param.OutGeometryType = GeometryType.Polygon;

            Polygon horizon = GeometryEngine.Instance.GeodesicEllipse(param,
                                                                      sr)
                                                                      as Polygon;

This is a sample of the failing spatial references:

EPSG:2044, Hanoi_1972_GK_Zone_18

EPSG:2065, S-JTSK_Ferro_Krovak

EPSG:2081, Chos_Malal_1914_Argentina_2
EPSG:2082, Pampa_del_Castillo_Argentina_2
EPSG:2083, Hito_XVIII_1963_Argentina_2

EPSG:2091, South_Yemen_GK_Zone_8

EPSG:2133, NZGD_2000_UTM_Zone_58S

EPSG:2177, ETRS_1989_Poland_CS2000_Zone_6

EPSG:2193, NZGD_2000_New_Zealand_Transverse_Mercator
EPSG:2195, NAD_1983_HARN_UTM_Zone_2S

EPSG:2206, ED_1950_3_Degree_GK_Zone_9

EPSG:2308, Batavia_TM_109_SE
EPSG:2309, WGS_1984_TM_116_SE

EPSG:2315, Campo_Inchauspe_UTM_19S

EPSG:2327, Xian_1980_GK_Zone_13

EPSG:2349, Xian_1980_3_Degree_GK_Zone_25

EPSG:2401, Beijing_1954_3_Degree_GK_Zone_25

EPSG:2523, Pulkovo_1942_3_Degree_GK_Zone_7

EPSG:2550, Samboja_UTM_Zone_50S

EPSG:2736, Tete_UTM_Zone_36S

EPSG:3318, IGC_1962_Congo_TM_Zone_12
EPSG:20248, AGD_1966_AMG_Zone_48

EPSG:22171, POSGAR_1998_Argentina_Zone_1

EPSG:22191, Argentina_Zone_1

EPSG:23877, DGN_1995_UTM_Zone_47S

EPSG:23886, Indonesian_1974_UTM_Zone_46S

EPSG:24877, PSAD_1956_UTM_Zone_17S

EPSG:28348, GDA_1994_MGA_Zone_48

EPSG:29177, SAD_1969_UTM_Zone_17S

EPSG:31266, MGI_3_Degree_Gauss_Zone_6

EPSG:32301, WGS_1972_UTM_Zone_1S

EPSG:32701, WGS_1984_UTM_Zone_1S

EPSG:32766, WGS_1984_TM_36_SE

Any suggestions appreciated.

[AnnetteLocke · Accepted solution]

Alan,

It might just be a copy/paste issue, but your formula for the circumference has a typo. 3 * (+b) should be 3 * (a + b). Another issue is that the (0, 0) may not be in the horizon of the projected coordinate system, so you can't use it as the center point if you are creating the geodesic ellipse in the projected coordinate system. However, you might be able use (0, 0) as the center point if you create the geodesic ellipse in the geographic coordinate system and then project it to the projected coordinate system. I tested this for a few of the spatial references that were failing, and it seems okay. 

double a = sr.Datum.SpheroidSemiMajorAxis;
double b = sr.Datum.SpheroidSemiMinorAxis;
double circumference = Math.PI * (3 * (a + b) - Math.Sqrt((3 * a + b) * (a + 3 * b)));

GeodesicEllipseParameter param = new GeodesicEllipseParameter();

param.LinearUnit = sr.Unit as LinearUnit;

param.SemiAxis1Length = circumference / sr.Unit.ConversionFactor;
param.SemiAxis2Length = circumference / sr.Unit.ConversionFactor;
param.AxisDirection = 0;
param.Center = new Coordinate2D(0, 0);

param.VertexCount = 1000;
param.OutGeometryType = GeometryType.Polygon;

SpatialReference gcs = sr.Gcs;
Polygon gcsPoly = GeometryEngine.Instance.GeodesicEllipse(param, gcs) as Polygon;
Polygon pcsPoly = GeometryEngine.Instance.Project(gcsPoly, sr) as Polygon;
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[AlanStewart]

Annette,

Thanks so much for looking at this. Yes that was a transcription error, I should have noticed it. I'm not a math or spatial reference specialist, I doubt I ever would have come up with this solution on my own. It does successfully handle the entire set of spatial references I have tested against.

Thanks again!