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

[Lesparre, Jochem]

Charles,

> My results use the *ellipsoidal* gnomonic projection.  

Ah, that explains the accurate results.


> This [spherical] is still 10 times better that your RD results.

That's why I suggested not to use point in polygon in RD with long segments.


> 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.  

Indeed, but most border polygons have no segments longer than 200 m (the Netherlands is a very crowded country). 


> 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.

It would indeed. This is tricky however. The digital storage system doesn't know whether the legal authorities meant to define a geodesic, rhumb line or a y=ax+b line in the national projected CRS (RD) for a policy or a zoning border. To avoid this problem, it is more safe to force the legal authorities to be specific, without having to explain them the difference between a geodesic, rhumb line, projections etc. 

Regards, Jochem


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
>>



Disclaimer:
De inhoud van dit bericht is uitsluitend bestemd voor geadresseerde.
Gebruik van de inhoud van dit bericht door anderen zonder toestemming van het Kadaster 
is onrechtmatig. Mocht dit bericht ten onrechte bij u terecht komen, dan verzoeken wij u 
dit direct te melden aan de verzender en het bericht te vernietigen. 
Aan de inhoud van dit bericht kunnen geen rechten worden ontleend.

Disclaimer:
The content of this message is meant to be received by the addressee only.
Use of the content of this message by anyone other than the addressee without the consent 
of the Kadaster is unlawful. If you have received this message, but are not the addressee, 
please contact the sender immediately and destroy the message.
No rights can be derived from the content of this message