ST_Distance, use_spheroid default

[gitai]

Per docs:

​http://postgis.net/docs/ST_Distance.html

"..For geography type defaults to return spheroidal minimum distance between two geographies in meters. ..For geography type defaults to return the minimum distance around WGS 84 spheroid between two geographies in meters. Pass in false to return answer in sphere instead of spheroid."

Test:

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326));

— output (looks like a bug): 0.234024137979163

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326), false);

— output: 19967.534156972

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326), true);

— output: 20013.836017253

Env:

Win 8.1

pgAdmin III v1.18.1

PostgreSQL 9.3.5, compiled by Visual C++ build 1600, 64-bit (EDB)

POSTGIS="2.1.3 r12547" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER

[nicklas]

What is happening is that ST_MakePoint won't give you a geography type but a geoemtry type ​http://postgis.net/docs/manual-dev/ST_MakePoint.html

Then I guess when you call ST_Distance with false/true for Sphere you will get an autocast to geography since that option is not valid for geometry type.

I guess the ticket as it stands should be closed as invalid, but maybe it is worth discussing if the auto cast in this case is a little dangerous.

[gitai]

hmm.. new to postgis .. very confusing..

​http://postgis.net/docs/manual-dev/ST_MakePoint.html

1. "Note x is longitude and y is latitude"

—Return point marked as WGS 84 long lat

2. SELECT ST_SetSRID(ST_MakePoint(-71.1043443253471, 42.3150676015829),4326);

Seems very dangerous .. Search the internet and you will see that ST_MakePoint is very popular method to quickly create points for params for geography calculations.. I guess people rely on the autocast.

BTW, Is there a way to create a geography in PostGIS without cast?

[nicklas]

​http://postgis.net/docs/manual-dev/ST_GeogFromText.html

[nicklas]

But the easiest in your case I guess would be to do the cast manually:

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326)::geography, ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326)::geography);

[gitai]

Actually, coming to think of it, the bug is that the ST_Distance interface does not properly enforce type safety between geometry and geography.

[gitai]

Since making points is such a common scenario, I think it would be better to have a distinct method (e.g ST_MakeGeogPoint).

[nicklas]

Replying to gitai:

Actually, coming to think of it, the bug is that the ST_Distance interface does not properly enforce type safety between geometry and geography.

Well, it is a result of function overload.

It behaves different depending on what type you are giving to it.

That is very common behavior, but what might be confusing is the distinction between geometry and geography type.

​http://postgis.net/docs/manual-dev/using_postgis_dbmanagement.html#PostGIS_Geography

A geometry type and geography type can both hold spatial information in srid 4326. But the calculations on geography type will be on the sphere and on geometry type will be planar.

I wouldn't call it a bug since it is by design, but maybe confusing first time you see it.

[gitai]

I'm new to postgis, not gis … and no your wrong, the function behaves differently as result of an overload *as well as* an implicit default parameter - which makes it a tricky interface.

Furthermore, your explanation does not seem consistent with how ST_DWithin works:

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326));

— output: 0.234024137979163

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326)::geography, ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326)::geography);

— output: 20013.836017253

select ST_DWithin(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), 20 * 1000)

— output: t

[gitai]

Also for non-edge case:

select ST_DWithin(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326), 20 * 1000)

— output: t

[gitai]

I expect to get an error passing a geometry to a geography interface (e.g use_spheroid overload), which would surface the mix up instantly.

[nicklas]

I don't get your point about ST_Dwithin. When the function uses planar calculations o lat lon of course the geometries is within 20000 units.

[gitai]

Actually, I misunderstood ST_DWithin to be at least… sorry about that. Your explanation is indeed consistent with it's definition.

But I still hold my argument for ST_Distance

[nicklas]

Replying to gitai:

I expect to get an error passing a geometry to a geography interface (e.g use_spheroid overload), which would surface the mix up instantly.

Yes, I see your point here. This is what I meant earlier about discussing if this autocast is confusing things more than doing good. (I actually just guess that it is an autocast that is happening from what you describe. Haven't had time to look at it).

Let's wait for the gurus to comment

[pramsey]

Removing the implicit geometry::geography cast would be a big change, though not probably one that would effect many people, since nearly every geography function has an exact geometry analogue, so the case of calling st_geogfunction(geom) almost never happens current. The reverse, on the other hand, happens quite a lot.

And, in fact, we don't have an implicit cast in the other direct, for reasons that are probably good, but I cannot recall.

select st_x('POINT(1 2)'::geography);
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

[pramsey]

That said "remove implicit geometry::geography type cast" should be filed as an enhancement, not a bug.

[gitai]

The ticket is about ST_Distance, not about removing the implicit geometry::geography cast across the board… fact is that calling st_geogfunction(geom) is exactly what (unintentionally) happened to me due to function overloading. Point being is that there is no "hint" as to what ST_MakePoint returns as opposed to functions like ST_GeogFromText where its clear, I mean like what makes a point more special for geometry than geography? Now I'm not suggesting to rename ST_MakePoint to ST_MakeGeomPoint,which is what it really does, as this would probably be a breaking change, but at least enforce the type safety in ST_Distance in a consistent manner.. I don't understand how this can be classified as "enhancement" seems like as a small overlooked usability bug that's painful to spot especially for newcomers to postgis (like me

[robe]

gitai — I don't think there is a way to satisfy this request short of what Pramsey is saying of removing the autocast in geometry::geography which many have complained about so would be a good enhancement in my book.

Both geometry and geography have ST_Distance functions (that take no additional args but the spatial objects) and geography is the only one to take the additional use_spheroid option which is how you managed to fall into the lucky trap.

So when you use the 2 argument one since you have 2 geometries it goes into the best fitting function (that ST_Distance that takes geometries). When you pass it the additional arg, only the geography version has the extra arg — and geometry says — hey I can pretend to be a geography to get myself in this door. By the time it gets into the function door, it already has been cast to a geography so there is no way we can provide a hint to say "hey dude you just morphed yourself into a geography".

So the fact you happened to waltz into the right geography function when passing your geometries in is nothing we can protect you from short of removing the autocast or renaming our functions to have Geom / Geog in name. Such a renaming would annoy a lot of people. However I do think many people (I included) would like the proposed get rid of the autocast. The autocast bytes you in very subtle ways like when when doing && between a geometry and geograpy

[gitai]

Replying to robe:

So the fact you happened to waltz into the right geography function when passing your geometries in is nothing we can protect you from short of removing the autocast …

Which is exactly the scope of the ticket … I'm not familiar with source to hold a discussion as to why removing the autocast from ST_Distance requires also removing autocast from all other functions, which appears to be why this issue is rejected, but probably there is reason..

[robe]

autocast is never done per function, its done with a cast object and automatically applies to all function. The thing that controls it is this:

CREATE CAST (geometry AS geography)
  WITH FUNCTION postgis.geography(geometry)
  AS IMPLICIT;

Note the term — IMPLICIT (that means that if there is a function for geography and no equivalent for geometry, then geometry can auto cast itself to geography to use it)

So if you want to experiment yourself with what pramsey is proposing do this:

-- this line needed if you installed postgis as an extension
ALTER EXTENSION postgis DROP CAST (geometry AS geography); 
DROP CAST (geometry AS geography);

 CREATE CAST (geometry AS geography)
  WITH FUNCTION geography(geometry)
  AS ASSIGNMENT;


--this line needed if you installed postgis as an extension
ALTER EXTENSION postgis ADD CAST (geometry AS geography); 

Once you do that, you'll see that this no longer works

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326), ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326), true); 

and you have to explicitly cast to

select ST_Distance(ST_SetSRID(ST_MakePoint(-73.9662500, 40.7834300), 4326)::geography, ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326)::geography, true); 

As a general rule, all constructor functions (except for the ones that have ST_Geog.. or that explicitly take in input as geography return geometry (not geography).

[gitai]

Replying to robe:

autocast is never done per function, its done with a cast object and automatically applies to all function. The thing that controls it is this:

CREATE CAST (geometry AS geography)
  WITH FUNCTION postgis.geography(geometry)
  AS IMPLICIT;

Ah.. Implicit type conversion.. That’s indeed a different story.

}}}

Note the term — IMPLICIT (that means that if there is a function for geography and no equivalent for geometry, then geometry can auto cast itself to geography to use it)

I assume this was implemented in order to maintain a single set of constructor functions; however, removing the autocast will probably result in more explicit casts.

I personally don’t think this is less annoying than having two sets of constructor functions one for geography and one for geometry.. You actually need to type less by convention plus it’s one less “additional concern” (e.g cast):

Geography construction (Geometry function call + cast)

ST_SetSRID(ST_MakePoint(-74.1960983, 40.739418), 4326)::geography

Vs

Geography construction (Geography function call)

ST_SetSRID(ST_MakeGeogPoint(-74.1960983, 40.739418), 4326)

BTW, just out of curiosity do constructer functions return geometries because geometry usage is more frequent than geography?

[robe]

geometry was the first and only before, some constructors are also simply much harder to implement in geography, and geometry is still the most popular. So yes.

Also while it may be easier for users, we didn't want to bother doing this for all geometry functions that could be used in this way because we already have too many functions for people to get confused with. It very much clutters the code base. An ST_MakePoint that does take as input an SRID would be nice though. Anyway this is turning into a discussion rather than a ticket. And strk the "No discussions in tickets" hound I expect will come barking soon.

Best to continue this on PostGIS development or PostGIS users list. You should join if you haven't already:

​http://postgis.net/support

I imagine a lot of people feel the way you do and it would be a good discussion topic.

[gitai]

Thanks!