[GRASS-dev] latlong problem

[Glynn Clements]

Martin Landa wrote:

> > > there is a significant problem in LatLong location. Try the following:
> > >
> > > g.region n=80 s=-80 w=-170 e=170 res=1;
> > > v.in.region out=reg;
> > > g.region n=90 s=-90 w=-180 e=180 res=1;
> > > d.erase;d.vect reg
> > > [1] should be [2]
> > >
> > > I prepared a simple patch to fix this bug.
> > >
> > > I am not sure whether it is OK. Anyway testing "(east1 - east2) > 180"
> > > does not make sense to me... or am I wrong?
> >
> > I've committed a change to v.in.region to add additional vertices if a
> > lat/lon region is wider than 179°, so such regions shouldn't get
> > turned inside-out by code which uses the "shortest route" maxim.
> >
> > Also, I've modified D_polygon_clip (used by "d.vect render=c") to
> > handle longitude wrap. It doesn't handle polygons which contain a
> > pole, but the other rendering methods don't handle that either (in
> > fact, I doubt that many things *do* handle that case). OTOH, it should
> > handle polygons which make multiple "orbits" (e.g. an "apple peel"
> > spiral).
> 
> g.region n=90 s=-90 w=-180 e=180;
> v.in.region out=reg
> 
> works well, thanks!
> 
> Trying to display world map, the problem with Antarctica still
> remains, now it is not filled. Maybe "buggy" boundaries. Not sure.
> Just for comment, when I enable (east1 - east2) > 360, the area is
> filled correctly.
> 
> In summary,
> 
> render=g -> Antarctica filled wrongly
> render=r -> OK (don't d.zoom -p)
> render=d -> as render=r
> render=c -> now not filled

The changes to handle longitude wrapping are likely to discard
polygons which bound a pole.

More precisely: the euclidify() function traverses the list of
longitude values; if a value differs from the previous one by more
than ±180°, it has the appropriate multiple of 360° added to it to
bring it into range, so that all line segments span no more than 180°. 

If the gap between the first and last values exceeds 180°, the
(presumed) polygon cannot be mapped to the Euclidean plane.

One way to look at is: if you adjust the first point to fall within
±180° of the final point, you then have to adjust the second point to
fall within ±180° of the first point, and so on, resulting in an
infinite loop.

Another way to look at it is that there is no way to cut the
cylindrical lat/lon surface so that it forms a plane without breaking
the polygon's boundary.

> The testing location is here if you are interested
> 
> http://gama.fsv.cvut.cz/~landa/.../grass/latlong.tar.gz

I note that the largest polygon for Antarctica contains a "fix-up"
segment between the vertices:

 179.9999     -89.9999    
 -179.9999    -89.9999    

If I suppress the "shortest path" maxim for segments where |lat|>89.9°
(which will typically be fix-up for the pole), the polygon is
considered closed.

I also needed to coerce the initial values into the range
-180°<lon<180° (rather than 0°<lon<360°) to avoid leaving a blank
triangle. Essentially, the radii which (along with the -89.9999th
paralell) formed the cut were both being placed at the left-hand edge. 

With those two changes, Antarctica rendered correctly. Hopefully,
that's likely to be good enough for most practical purposes.

-- 
Glynn Clements <glynn at gclements.plus.com>