[PROJ] How does proj deal with ellipsoid with respect to reprojection

[Charles Karney]

Jochem,

My results use the *ellipsoidal* gnomonic projection.  (This isn't
provided by PROJ.  GeographicLib does provide it.  For my experiments I
used the MATLAB/Octave version of GeographicLib.)  With the spherical
version of the gnomonic projections, the deviations are about 170 times
larger (deviation is 166.3 m for segment length 500km and it varies
quadratically with the length of the segment).  This is still 10 times
better that your RD results.

Requiring that all polygon segments be 200 m or less would seem incur a
large recurring cost for storing and processing (and even understanding)
the data.  Surely a better scheme would be to store the data with the
minimum number of vertices but (optionally) computing intermediate
vertices on the fly when delivering the data to the end user.

   --Charles

On 3/30/20 4:12 AM, Lesparre, Jochem wrote:
> Charles,
> 
> Thanks, that is very interesting! It shows that the gnomic projection gives quite good results, despite neglecting the flattening. For point in polygon with 1 mm accuracy, segments of 50 km are allowed. This makes it probably more efficient than adding additional points every 200 m for +proj=lonlat. Since geodetic knowledge is becoming more and more scarce in the Netherlands, I probably keep advising the use of segments of maximum 200 m to avoid misinterpretation by laymen.
> 
> Regards, Jochem
> 
> -----Original Message-----
> From: Charles Karney <charles.karney at gmail.com> On Behalf Of Charles Karney
> Sent: zondag 29 maart 2020 16:17
> To: Lesparre, Jochem <Jochem.Lesparre at kadaster.nl>; proj at lists.osgeo.org
> Cc: Pierre Abbat <phma at bezitopo.org>
> Subject: Re: [PROJ] How does proj deal with ellipsoid with respect to reprojection
> 
> Jochem,
> 
> Here's the equivalent table for the ellipsoidal gnomonic projection as it applies (roughly) to the Netherlands.  I use the ellipsoidal gnomonic projection centered at lat = 52 deg.  I considered geodesic line segments of length S whose centers are 500 km from the center of the projection.  I then computed D, the maximum deviation of the straight line in the projection from the geodesic.
> 
>         S     max-deviation
>       500 km   99 mm
>       200 km   16 mm
>       100 km    4 mm
>        50 km  985 um
>        20 km  158 um
>        10 km   40 um
>         5 km   10 um
>         2 km    2 um
>         1 km  394 nm
>        500 m   99 nm
>        200 m   18 nm
> 
> On 3/29/20 2:28 AM, Lesparre, Jochem wrote:
>> Doing point in polygon in a projection will result in occasional wrong
>> conclusions. A point near de edge can seem to be inside the polygon
>> while it's outside, or the other way around, since a straight line in
>> the projection deviates from the geodesic.
>>
>> I analysed this problem for the Netherlands (51 - 55 degrees north) in
>> the azimuthal projection (+proj=sterea) of the national coordinate
>> reference system RD (epsg:28992) and in plate-caree projection
>> (+proj=lonlat). The deviation depends on the length, orientation and
>> location of a polygon segment. I computed the maximum possible
>> deviation in the Netherlands for both projections to advise the Dutch
>> government not to allow any segments longer than 200 m in a new
>> digital storage system for policy and zoning borders.
>>
>> Table for the Netherlands:
>> Segment length, RD deviation, lonlat deviation;
>>
>> 500 km, 160 m, 6 km;
>> 200 km, 25 m, 1 km;
>> 100 km, 6.4 m, 0.2 km;
>> 50 km, 1.6 m, 60 m;
>> 20 km, 26 cm, 9.7 m;
>> 10 km, 8 cm, 2.4 m;
>> 5 km, 3 cm, 60 cm;
>> 2 km, 5 mm, 9.7 cm;
>> 1 km, 1.3 mm, 2.4 cm;
>> 500 m, 0.3 mm, 6 mm;
>> 200 m, <0.1 mm, 1 mm;
>> 100 m, <0.1 mm, 0.2 mm;
>> 50 m, <0.1 mm, <0.1 mm
>>
>> This means that when you first split long segments of polygons by
>> adding enough points along the geodesic, you could do a point in
>> polygon even in lonlat. The required distance between the added points
>> depends on the location (latitude). Near the poles it will get a bit
>> tricky, as the maximum deviation increases to 50% of the size of the segment length.
>>
>> Regards, Jochem
>>