When a SRID is provided, use its units for metrics of length, distance, and area instead of input units

[roji]

Does NetTopologySuite include any sort of support for geography (as opposed to geometry)? I'm specifically interested in communicating with PostgreSQL/PostGIS.

If not, is geography support something that is planned or possible? Does JTS implement it?

[FObermaier]

NTS does not perform any special geography operations.
The workflow for handling geography would be to

• project the geography to an adequate projected spatial reference system,
• do the operation/predicate wanted and
• reproject the result back to geography (WGS84)

JTS doesn't implement specific geography functions or operations either.

[roji]

Thanks for your answer.

For supporting PostGIS's geometry and geography data types (which are two distinct types), no actual special support is needed from NTS since both types have exactly the same representation in both WKT and WKB. So NetTopologySuite's Point can be used to represent both a geometry point or a geography point.

[airbreather]

I've started a bit on this: https://github.com/NetTopologySuite/NetTopologySuite.Geography

I just sorta went through some of this page for SQL Server and added stuff (not everything yet, but all OGC instance methods are represented).

Where there were direct equivalents on GeoAPI.Geometries.IGeometry, I (tried to) make the signature the same, otherwise I used my own judgment.

[airbreather]

Packages are now being auto-built and auto-pushed to my personal MyGet feed for now: https://www.myget.org/feed/airbreather/package/nuget/NetTopologySuite.Geography

[airbreather]

I'm unassigning myself from this one and waiting for the EF Core team to gather community feedback and provide a recommendation for how it makes sense to move forward with geography from their perspective.

Not that I'm saying that EF Core is the authoritative source for this, but I expect them to be a very large use case, and based on other communication, there may very well be more appropriate alternative ways to support this kind of thing client-side (even if it makes EF Core's SQL Server database provider a pain).

[bricelam]

Sorry! I meant to follow up on this yesterday. I just merged dotnet/efcore#13408 which allows you to map to SQL Server geography columns by specifying the column type:

class SpatialThing
{
    public int Id { get; set; }

    [Column(TypeName = "geography")]
    public Point Location { get; set; }
}

This will alter the translation of a few members: (note, the semantics are similar to how WKT is handled)

NTS	Transact-SQL
X	Long
Y	Lat
ExteriorRing	RingN(1)
NumInteriorRings	NumRings - 1
GetInteriorRingN(n)	RingN(n + 2)

The members that are't available in SQL Server will be evaluated on the client using the NTS implementations. This will probably be unexpected, but if you know you're working with geography, you probably won't call these methods anyway.

We're considering making geography the default mapping since this seems to be what most people want, and you't have to explicitly map to geometry.

Again, all this is open to change based on feedback.

[airbreather]

The decisions that EF Core made on this issue for 2.2 are the best I could also come up with:

1. geography is the default mapping for spatial types
   • I definitely did not expect to support this idea when we first started talking about this, since geometry more directly maps to the concepts that NTS / GeoAPI.NET represent, but upon reviewing everything, I came to realize that there's actually very little difference, and geography is likely to be significantly more useful for real-world GIS applications.
   • That gap is actually just representative of a broader observation that I have about the state of existing GIS software packages, which I unproductively rant about later in this comment.
2. EF Core's SQL Server adapter handles bridging NTS's and SQL Server's incompatible representations of compatible concepts, with some help from an I/O bridge which we now maintain.
3. Client-side computation will use NTS, so some length / area / etc. measurements will not be accurate without something changing on our part.

That third bullet point is why I still want to keep this issue open.

As annoying as it is to suggest major deviations from JTS, I think it's inappropriate for us to keep this situation working the way it does:

IPolygon poly = GetGeometry();
Console.WriteLine(poly.SRID); // writes 4326
Console.WriteLine(poly.Area); // writes a useless value

This is especially true now that our interfaces are used to enable this kind of scenario where that's exactly how to get a useful area (sorry for what's probably bad hygiene in this LINQ query):

var q = from q1 in GetQueryableFeatures1()
        from q2 in GetQueryableFeatures2()
        where q1.Geom.Intersects(q2.Geom)
        select (id1: q1.Id, id2: q2.Id, area: q1.Geom.Intersection(q2.Geom).Area);
foreach (var (id1, id2, area) in q)
{
    Console.WriteLine($"{id1} and {id2} share {area} square meters.");
}

I think that "real-world metric" operations in NTS (distance / area / length / more?) should start using SRID to tell us how to compute it. Ideally, using ellipsoid math instead of projecting, but the realist in me is guessing that it's going to use projections for the first round, and then there won't be another round.

Unproductive Rant

To elaborate on that issue I have with the state of existing GIS software packages: I know of precious few software packages (especially FOSS) that really "understand" that we're dealing with an ellipsoidal representation of the Earth where a point is identified by a pair of angles.

• For the most part, packages just pretend like the lat/lon values are planar coordinates because that math is easier (therefore faster) and equally correct in many useful cases.
• Whenever that charade falls apart, the go-to solution seems to be to rewrite all the individual coordinates using a projection so that the planar math works, and hope that nothing crosses a particular line of longitude (exactly which line of longitude may depend on how the projection is defined).
  • So if we're using an equal-area projection, we also hope that the input has a high enough resolution to minimize the effects of distortion, but not so high of a resolution that minor inaccuracies resulting from using finite-precision arithmetic make strange things happen.
• Even dealing with regular ol' unprojected coordinates in ways that are usually correct, we still run into problems inherent to the disconnect between our model and reality. For example, the Envelope of a MULTIPOLYGON that represents the U.S. state of Alaska also covers Moscow, Dublin, all of Estonia, all of Ireland, etc. Careless otherwise-valid data processing could make this worse than just the Envelope.
  • The U.S. is probably not going to get Alaska to give part of the Aleutian islands to Russia, and I doubt that anybody's going to erect an impenetrable force field all along the antimeridian stopping anybody from crossing it (the U.S. is currently having a hard enough time building a comparatively insignificant barrier along part of a land border), so it feels like everybody is just hoping that nobody really cares about any events that cross this imaginary line.

So, forget about accurately modeling (non-geostationary) satellites that orbit the Earth constantly in one direction. For all we care, they actually follow Pac-Man logic: once they reach a particular point in their orbit, they just zap over to "the other side" to continue their journey.

If I want a geometry that represents the water area within 1 mile of the coast of any land on Earth, then it's not just LAND.Buffer(1 mile).Difference(LAND); I'm pretty sure that it's actually:

• Split LAND into components that our projection will keep contiguous
  • Probably two: one for everything between 45-135 longitude and one for the rest
• Do the buffer(1 mile) on each component
• Unproject each component using the original projection
• Cut any individual components that now happen to cross the antimeridian
• LAND_AND_SOME_WATER = Union the results
• Your result is LAND_AND_SOME_WATER.Difference(LAND)

Of course, what can be done about this? I'm sure that I'm not the first, nor will I be the last, to be bothered by all this, but it's not like we got here by incompetence or anything. It's just really hard to model that, and what we have is "good enough" for enough real-world use cases that developers just live with the pain and extra grey hairs of working around the flaws of imperfect solutions.

[roji]

Loved the rant :)

Just one point - geography is the default for SQL Server, but I don't think this constitutes a standard EF Core decision. For example, Npgsql maps to geometry by default, because PostGIS geography is significantly more limited and so a bit useless, depending on exactly what you're trying to accomplish. It's more accurate to say that the mapping is completely database-dependent.

[airbreather]

Ideally, the magic should be opt-out, which, for the first round, pretty much means that NTS needs to be able to access ProjNet stuff.

@FObermaier has already expressed in NetTopologySuite/ProjNet4GeoAPI#44 (comment) that ProjNet should not reference anything from NetTopologySuite, which is perfect, since that lets NetTopologySuite reference ProjNet 😀.

So the basic plan of attack to get to a place where it's definitely more useful than what we do today:

1. Move the logic for at least the SRID-dependent operations into something that can be injected, initially injecting an implementation that works exactly how Geometry works today.
2. Remove ProjNet --> NetTopologySuite and ProjNet --> GeoAPI references.
3. Add NetTopologySuite --> ProjNet reference
4. Build an implementation of those SRID-dependent operations that implements them as:
   • Passthrough when the SRID identifies a projected coordinate system,
   • Run it on a copy that's been transformed to some projected coordinate system otherwise.

Bad things about that:

1. Geometry instances can be really heavyweight... copying + coordinate transformation alone is probably going to take the bulk of the time in the bulk of the cases where performance matters.
2. On the flip side, if the geometry instance has "too few" points, then I think we'll accumulate maybe "too much" error by calculating based on perfectly straight lines between vertices?
   • I haven't ever tested this hypothesis rigorously, so it's entirely possible that pointwise transformation is perfectly fine. If this particular hypothesis I have is correct, then it would be possible to see small* differences in results by re-testing after adding collinear points between unprojected vertices.
   • *these differences would be significantly more than the trivial differences we'd expect to result from imprecise floating-point calculations.
3. Vertices of projected "holes" might theoretically move outside of the polygon's projected "shell", even for perfectly valid input polygons, and I haven't been able to prove that it's impossible for area to go negative in cases where this happens.

Resolving those "bad things" would probably require us to implement the calculations using ellipsoidal math, using ProjNet for nothing more than its ability to extract the parameters we need for a given SRID, which seems very time-consuming.

[airbreather]

Resolving those "bad things" would probably require us to implement the calculations using ellipsoidal math, using ProjNet for nothing more than its ability to extract the parameters we need for a given SRID, which seems very time-consuming.

Then again, if all of it just flows from distance and area...

Distance can be calculated using Vincenty's formulae (Mike Gavaghan did a C# implementation which I've already copied to GitHub with permission).

https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid#Area_of_a_geodesic_polygon doesn't look too awful.

[airbreather]

Distance can be calculated using Vincenty's formulae

https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid#Area_of_a_geodesic_polygon doesn't look too awful.

Looks like the way to go is actually Charles F. F. Karney's Algorithms for geodesics + addenda. There's already an implementation, so it should just be a matter of porting it and/or updating suryapratap/GeographicLib or Flitesys/GeographicLib or oldrev/GeographicLib or similar.

[airbreather]

So, since we're probably going to wind up doing this via ProjNET's pointwise stuff, I've gone ahead and started a branch wip/sequence-transform that adds a NTS --> ProjNET reference and adds some stuff to Geometry and CoordinateSequence that makes it work. Not on develop yet because I think the ProjNET API could use a bit of a tune-up as well...

Notably, PackedDoubleCoordinateSequence and DotSpatialAffineCoordinateSequence get to use the Span<T> methods on MathTransform, because their coordinate data is laid out appropriately.

[airbreather]

2.0 is releasing without this, but I would like a future 2.x release to add support for this, on an opt-in basis to avoid breaking existing consumers.

[juliusfriedman]

Before creating a new issue I wanted to chime in here because this seems related.

I also read This But didn't find it helpful because even after projecting or using LengthLocationMap I get the same distance 0.00906611787394443

Given a LineString:

LINESTRING (-1.0858652591705322 52.70209884643555, -1.0853074789047241 52.701202392578125, -1.0851855278015137 52.70097732543945, -1.0850050449371338 52.70076370239258, -1.0846827030181885 52.70051193237305, -1.0843605995178223 52.70035171508789, -1.0840034484863281 52.700233459472656, -1.0836650133132935 52.700164794921875, -1.083327054977417 52.700130462646484, -1.0829579830169678 52.700138092041016, -1.0826040506362915 52.70017623901367, -1.082174301147461 52.70027160644531, -1.0816984176635742 52.700477600097656, -1.0813709497451782 52.70071792602539, -1.081089735031128 52.700984954833984, -1.0808309316635132 52.70120620727539, -1.0802907943725586 52.701622009277344, -1.0797139406204224 52.70204544067383, -1.0785545110702515 52.702880859375)

The Length = 0.00906611787394443
The StartPoint.Distance(EndPoint) = 0.0073524541496857278

SQL Server gives a length of 759.743066435061

How can I get to that value without having to manually select it from SQL Server using FromSql? Is there a way to convert easily?

[juliusfriedman]

I tried the following:

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

            var transformationFactory = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
            var transform = transformationFactory.CreateFromCoordinateSystems(
                originCoordinateSystem, targetCoordinateSystem);

            var pointA = ss.StartPoint;
            var pointB = ss.EndPoint;

            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}");

But the distance is 603.137021134185 which is still off of what SQL Server Reports for the given line string

if I modify slightly with the code in NetTopologySuite/ProjNet4GeoAPI#36
e.g. using the CalcUtmZone Then I get 501.69676329929422.

So in short I am totally lost here

[FObermaier]

You need to transform the whole geometry/linestring and not just start- and end-point.
As-is you are computing the point to point distance between start- and end-point

[juliusfriedman]

I see....

With this I get closer:

 var originCoordinateSystem = GeographicCoordinateSystem.WGS84;
            var targetCoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(CalcUtmZone(ss.StartPoint.X), true);

            var transformationFactory = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
            var transform = transformationFactory.CreateFromCoordinateSystems(
                originCoordinateSystem, targetCoordinateSystem);

            CoordinateList list = new CoordinateList();

            double distanceSum = 0;

            GeoAPI.Geometries.Coordinate lastPoint = null;

            for(int i = 0, e = ss.NumPoints; i < e; ++i)
            {
                var oldCoord = ss.GetCoordinateN(i);
                var newPoint = transform.MathTransform.Transform(new GeoAPI.Geometries.Coordinate(oldCoord.X, oldCoord.Y));

                if (lastPoint != null)
                {
                    distanceSum += newPoint.Distance(lastPoint);
                }

                lastPoint = newPoint;

            }

I get 759.59573359987223 as the distanceSum which is not much different than 759.743066435061 as returned by sql server.

Thank you!

BTW, looking forward to seeing this type of functionality be enhanced in future versions, I already specify the SRID and I feel it would be beneficial if there was at least a built in function set to tackle these types of problems especially when you are dealing with many LineStrings and many Points as you have to loop over every coordinate, transform and then calculate which can eat up memory vs having that information already provided via the SRID and easily available at the server level.