Distance between points with transform

[eiredrake]

Unfortunately I'm getting a little frustrated. I am using .net standard 2.0 with GeoAPI 1.7.5, NettopologySuite 1.15.1 and ProjNet4GeoAPI 1.4.1

I am trying to get the distance between two latitude/longitude points in meters that I know to be roughly 10km apart in the united states. All of my calculations are done using the same two GeoAPI.Geometry.IPoint classes with correctly set X/Y values.

I have done the calculation with postGIS in a query window and gotten: 9052.83580248723 meters (using SRID 4326 to 2163 ST_Transform() and ST_Distance())

I've done the calculation by converting using the Distance() method between the two points (i get: 0.089028956413068211 which is degrees if I understand correctly) and apply generic multiplier from radians distance of 111.325 wherein I get: 9911.1485726848186 m

I've done the haversine algo wherein I get 1449.7855112772756 m

I've also been sitting here for about an hour with google trying to get the distance via point transform using Proj4Net between 4326 wgs84 and pretty much anything that the WKT specifies a unit of meter and I'm getting wildly different answers none of which match. The code I'm using (along with the last WKT I tried) is below. All I want is the distance between two points in meters....I'm of the mind that using the tools I am using this shouldn't be all that complicated so clearly I must be doing something wrong. Can anyone kick me in the right direction?

var originWKT = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
var targetWKT = "PROJCS[\"US National Atlas Equal Area\",GEOGCS[\"Unspecified datum based upon the Clarke 1866 Authalic Sphere\",DATUM[\"Not_specified_based_on_Clarke_1866_Authalic_Sphere\",SPHEROID[\"Clarke 1866 Authalic Sphere\",6370997,0,AUTHORITY[\"EPSG\",\"7052\"]],AUTHORITY[\"EPSG\",\"6052\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4052\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",45],PARAMETER[\"longitude_of_center\",-100],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"2163\"]]";

var coordinateSystemFactory = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
var transformationFactory = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();

var originCoordinateSystem =  coordinateSystemFactory.CreateFromWkt(originWKT);
var targetCoordinateSystem = coordinateSystemFactory.CreateFromWkt(targetWKT);

var transform = transformationFactory.CreateFromCoordinateSystems(originCoordinateSystem, targetCoordinateSystem);

var pointACoordinate = new GeoAPI.Geometries.Coordinate(pointA.X, pointA.Y);
var pointBCoordinate = new GeoAPI.Geometries.Coordinate(pointB.X, pointB.Y);

var newPointA = transform.MathTransform.Transform(pointACoordinate);
var newPointB = transform.MathTransform.Transform(pointBCoordinate);

result = newPointB.Distance(newPointA);

[FObermaier]

Could you post the actual coordinates of pointA and pointB?

[eiredrake]

I can't show the original two points but I chose two arbitrary points roughly ten km from each other. I'm getting roughly the same result.

point a: POINT (-77.220888 39.168931)
point b: POINT (-77.150249 39.219669)

degreesDistance: 0.086972483953261043 (this I get just using the .Distance() method on the IPoint)
degreesToMeters: 9682.2117760967849
transform: 10701.365215911384 (wgs84 to Popular Visualisation CRS / Mercator)
haversine: 961.83508919368637 (i suspect there's something wrong with the formula here)
PostGIS: 8376.15057906345
meridianoutpost.com: 8.3 km
.movable-type.co.uk/scripts/latlong.html: 8.3km (haversine calculator)

[FObermaier]

I'm afraid Lambert_Azimuthal_Equal_Area projection is not supported (in targetWkt).
How did you get a result in the first place, I get NotSupportedException?

[eiredrake]

Oh I did. It's just one of the many projections I tried. I was going through them one by one copy/pasting the WKT of any that had a unit of 'meters' trying to get something that made sense.

[FObermaier]

var originCoordinateSystem = GeographicCoordinateSystem.WGS84;
var targetCoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(18, true);

var transform = transformationFactory.CreateFromCoordinateSystems(
    originCoordinateSystem, targetCoordinateSystem);

var pointA = new GeoAPI.Geometries.Coordinate(-77.220888, 39.168931);
var pointB = new GeoAPI.Geometries.Coordinate(-77.150249, 39.219669);

var pointACoordinate = new GeoAPI.Geometries.Coordinate(pointA.X, pointA.Y);
var pointBCoordinate = new GeoAPI.Geometries.Coordinate(pointB.X, pointB.Y);

var newPointA = transform.MathTransform.Transform(pointACoordinate);
var newPointB = transform.MathTransform.Transform(pointBCoordinate);

double result = newPointB.Distance(newPointA);
Console.WriteLine($"Distance between {pointA} and {pointB} = {result}");

produces

Distance between (-77.220888, 39.168931, NaN) and (-77.150249, 39.219669, NaN) = 8305,08674556441

[eiredrake]

Hrm....Getting a 'No support for transforming between the two specified coordinate systems' exception on the line where we call .CreateFromCoordinateSystems().

[eiredrake]

Also hrm... looking at what Universal Transverse Mercator coordinate system for zone 18.... (first time I've heard of it but google provides)

How do you determine what zone the points you're looking at are in? What if one point is in one zone and the second point is in another zone?

[FObermaier]

You can calculate the UTM zone using the following formula:

public static int CalcUtmZone(double lon) => (int) ((lon + 180.0)/6.0 + 1.0);

You need to make sure that the used projection is suitable for the input points.

SharpMap has a geospatial math utility class that might help

[eiredrake]

I have a similar class that calculates distance using similar methods. I'm just trying to understand how to do it using transforms. I had guessed that you transform both points into a projection that has it's units in meters and then just get the distance between them. Is that not how it works?

I wanna say I appreciate your help. I hate not understanding things.

[FObermaier]

To compute the distance between two projected coordinates, you need to make sure that

• both are projected into the same spatial coordinate system
• the used spatial coordinate system is suitable for the distance calculations (it does not make sense to use EPSG:31467 when your coordinates are in Australia)

[eiredrake]

right oranges to oranges.... but they are, aren't they? They're converted to WGS84 UMT... but then there's zone. Do you just put them in the same zone? or does that not matter because it's all WGS84 UMT?

[FObermaier]

A different UTM Zone is a different spatial reference system. The error might be negligable when you have neighboring UTM zones, but as a rule of thumb, avoid that.

See https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Overlapping_grids