Converting Geometry to a Geography yields an Invalid shape

311 Views Asked by At

I have two Geometry shapes of types Polygon or MultiPolygon [with Lat & Lng points] that I want to get their intersected shape. I do this as follows:

land_shape.STIntersection(grid_shape) as shape

Then, I check to see when

shape.STIsValid() = 1 and geography::STGeomFromText(shape.STAsText(), 4326).STIsValid() <> 1

and I get many results back. There are instances where the intersected shape as a geometry is valid, but as a geography isn't. This is an issue because I build on the geography shape to get the area.

I tested to see if converting both inputs (land_shape & grid_shape) to geography would also be invalid, but they are not invalid. So, inputs are not corrupt.

I'm not sure what could cause the intersected shape to be invalid as a geography. I tried using IsValidDetailed() and the reasons were:

24409: Not valid because some portion of polygon ring (1) lies in the interior of a polygon.
24413: Not valid because of two overlapping edges in curve (1).

But, why were they valid as a geometry then? I know that MakeValid() can fix it, but I just need to understand why that behavior happens to begin with, and if I can eliminate this issue without using SQL Server's MakeValid() as I read it isn't entirely accurate.

Update: Even after applying MakeValid() I still get very large areas returned How I calculate the area is as follows:

geography :: STGeomFromText(land_shape.STAsText(), 4326).STIntersection( geography :: STGeomFromText(grid_shape.STAsText(), 4326) ).STArea()

I thought initially that the problem was the Geography shapes being invalid, but looks like it might not just be that.

I tried applying ReorientObject() but still same issue

geography :: STGeomFromText(land_shape.STAsText(), 4326).ReorientObject().STIntersection( geography :: STGeomFromText(grid_shape.STAsText(), 4326).ReorientObject() ).STArea()
1

There are 1 best solutions below

8
Michael Entin On

Saving Geometry as WKT and reading it as Geography does not properly converts from geometry to geography. The two types have different interpretations of that WKT, even though they both use similar WKT and they use same coordinates. Edges between two vertices in Geometry are the straight lines on the flat map. Edges between the same two vertices in Geography follows the geodesic line - so it is a very different line.

Consider e.g. two loops:

  'linestring(30 40, 80 40, 80 60, 30 60, 30 40)'

and

  'linestring(50 42, 60 42, 60 50, 50 50, 50 42)'

On a flat map the second loop is strictly inside the first one - so you can make a polygon with first loop as shell, and second loop as hole - and it will be a valid polygon. But if you use the same WKT and try to interpret it as Geography - the second loop intersects the first one, and if you construct a polygon the same way - it is no longer valid.

enter image description here

Thus a polygon

 'polygon((30 40, 80 40, 80 60, 30 60, 30 40),(50 42, 60 42, 60 50, 50 50, 50 42))'

is valid as Geometry, but not valid as Geography. There are also issues around anti-meridian and poles. The proper conversion means approximating "flat" edges of Geometry with geodesic edges of Geography and vice versa - usually involving densification of the edge, i.e. adding some points in the middle. E.g. this is what Google BigQuery does in this case:

select st_geogfromtext(
    'polygon((30 40, 80 40, 80 60, 30 60, 30 40),(50 42, 60 42, 60 50, 50 50, 50 42))', 
    planar => TRUE)

note the planar argument, telling BigQuery the input comes from a flat geometry. The result is pretty large, but it follows the original "flat" shape with about 10 meter precision:

POLYGON((30 40, 30.1953125 40, 30.390625 40, 30.5859375 40, 30.78125 40, 30.9765625 40, 31.171875 40, 31.3671875 40, 31.5625 40, 31.7578125 40, 31.953125 40, 32.1484375 40, 32.34375 40, 32.5390625 40, 32.734375 40, 32.9296875 40, 33.125 40, 33.3203125 40, ....

One more issue is ring orientation - which does not matter for Geometry, but is very important for Geography, see https://learn.microsoft.com/en-us/sql/relational-databases/spatial/polygon?view=sql-server-ver16#orientation-of-spatial-data