[GRASS-dev] latlong problem
glynn at gclements.plus.com
Fri Mar 9 13:35:54 CET 2007
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
> > >  should be 
> > >
> > > 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
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
I note that the largest polygon for Antarctica contains a "fix-up"
segment between the vertices:
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
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>
More information about the grass-dev