#1527: vector projection over wrapping boundary is split
----------------------+-----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.0
Component: Vector | Version: svn-trunk
Keywords: | Platform: Linux
Cpu: x86-64 |
----------------------+-----------------------------------------------------
When I project a vector from a location thet overlaps with a wrapping
boundary of a map, the vector is split in two. For example, I have a
longitude latitude location covering the whole world, that wraps somewhere
in the Pacific (that is points on the eastern boundary are also on the
western boundary). When I project a vector constructed by doing a grid on
a Lambert equal area location centered somewhere in the Pacific, the grid
cells (that are, in that case, no more rectangular) that fall on the
boundary are cut in two, with warnings like:
{{{
WARNING: Number of centroids exceeds number of areas: 900 > 847
WARNING: Number of incorrect boundaries: 507
WARNING: Number of centroids outside area: 53
}}}
I attach a tarball to reproduce the issue (there is a png showing the
issue in the tarball too). The script to do everything (setup the lambert
location, and then the world map and do the png) is ./bug_wrap.sh.
#1527: vector projection over wrapping boundary is split
----------------------+-----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.0
Component: Vector | Version: svn-trunk
Keywords: | Platform: Linux
Cpu: x86-64 |
----------------------+-----------------------------------------------------
Comment(by mmetz):
Please try attached patch for v.proj. The bug is in the proj4 library
which automatically wraps to -180 - 180 if the output location is latlon.
This does not work with topological vectors and causes the observed error
messages. The attached patch tries to automatically figure out if wrapping
coordinates to the range 0-360 is needed and wraps if needed (should be
made a user option I guess). BTW, this bug appears also in GRASS 6.4, only
that 6.4 does less thorough topological error checking.
{{{
main.c:54: attention : ‘src_box.N’ may be used uninitialized in this
function
main.c:54: attention : ‘src_box.S’ may be used uninitialized in this
function
main.c:54: attention : ‘src_box.E’ may be used uninitialized in this
function
main.c:54: attention : ‘src_box.W’ may be used uninitialized in this
function
main.c:54: attention : ‘tgt_box.W’ may be used uninitialized in this
function
}}}
The bug seems to be fixed for my reproducer case. I will reopen if it
doesn't fix my 'real world' case.
Reopening the ticket because nothing is fixed yet. I would suggest to add
a new option for wrapping to 0,360 instead of -180,180, and maybe add a
test if wrapping would be needed. The purpose of the attached patch was to
demonstrate that the problem could be solved by wrapping to 0,360. But I
don't think this should be done automatically.
I have tested with the grass64_release svn branch and indeed the bug is
also present here. It is less visible, because the centroids are there
and the boundaries too, but centroids are not associated to the
boundaries. There is no error message but it can be seen that centroids
are not associated to areas:
{{{
Number of centroids: 900
Number of areas: 840
Number of isles: 7
Number of incorrect boundaries: 456
Number of centroids outside area: 60
}}}
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by pertusus):
When I look at the resulting map in bug_wrap/PERMANENT with the default
(world) region, there seem to be no centroid for everything coming from
the western boundary. I attach a tarball for 7.0 with this additional
image and the corresponding script.
Not sure if it was introduced by the patch or not.
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by pertusus):
A select doesn't cross the boundary. Only the eastern part is selected.
I do a grid on the world map and select it by the Lambert map. I have
attached a tarball bug_grass_proj_wrap_select.tar.gz to reproduce. The
difference between the Lambert map and the selected world map grid may be
seen by comparing world.png (the Lambert map) and selected.png (the world
map grid selected).
The bug_grass_proj_wrap_select.tar.gz also shows the issue mentioned above
on the missing visible centroids in the world image. I guess that both
issues are related.
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by mmetz):
Please try trunk r50023. There is a new option to disable longitude
wrapping directly in the proj4 library
Disabling longitude wrapping could (should?) be the default since usually
there is no difference, and if errors occur, the default wrapping needs to
be disabled in proj4 in order to preserve vector topology.
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by mmetz):
Replying to [comment:7 pertusus]:
> A select doesn't cross the boundary. Only the eastern part is selected.
>
> I do a grid on the world map and select it by the Lambert map. I have
attached a tarball bug_grass_proj_wrap_select.tar.gz to reproduce. The
difference between the Lambert map and the selected world map grid may be
seen by comparing world.png (the Lambert map) and selected.png (the world
map grid selected).
Not sure how to handle that. The result of v.select is ok if the grid is
created on a region with n=90 s=-90 e=360 w=0.
This would become a problem with e.g. world boundaries wrapped to -180,180
instead of 0,360. In this case the world boundaries would need to be
transformed with x shift = 360, the original and transformed world
boundaries patched and cleaned (preferably also clipped to 0,360) and then
used as input for v.select together with the Lambert map.
>
> The bug_grass_proj_wrap_select.tar.gz also shows the issue mentioned
above on the missing visible centroids in the world image. I guess that
both issues are related.
I don't think so. The missing centroids are caused by the display
routines.
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by pertusus):
Replying to [comment:10 mmetz]:
> Not sure how to handle that. The result of v.select is ok if the grid is
created on a region with n=90 s=-90 e=360 w=0.
But then, unless I have misunderstood something, the same issue will arise
for the projection of another vector that crosses the 0/360 border? Maybe
this should be discussed in another bug, but when one has vectors all over
the world and one wants to work with those vectors, it is not possible to
find a region in which the vectors do not overlap with a boundary. But I
may have missed something.
> This would become a problem with e.g. world boundaries wrapped to
-180,180 instead of 0,360. In this case the world boundaries would need to
be transformed with x shift = 360, the original and transformed world
boundaries patched and cleaned (preferably also clipped to 0,360) and then
used as input for v.select together with the Lambert map.
Unless I missed something this will work if the vector boundary does not
cross the 0/360 boundary, but will it work in a general case?
> > The bug_grass_proj_wrap_select.tar.gz also shows the issue mentioned
above on the missing visible centroids in the world image. I guess that
both issues are related.
>
> I don't think so. The missing centroids are caused by the display
routines.
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by mmetz):
Replying to [comment:11 pertusus]:
> Replying to [comment:10 mmetz]:
>
> > Not sure how to handle that. The result of v.select is ok if the grid
is created on a region with n=90 s=-90 e=360 w=0.
>
> But then, unless I have misunderstood something, the same issue will
arise for the projection of another vector that crosses the 0/360 border?
Maybe this should be discussed in another bug, but when one has vectors
all over the world and one wants to work with those vectors, it is not
possible to find a region in which the vectors do not overlap with a
boundary. But I may have missed something.
When projecting another vector that crosses the 0/360 border, you would
need again the new -w flag for v.proj and get a result similar to what you
have now, i.e. in the range 0,360.
Usually there are only 2 special cases when reprojecting vectors from
another CRS: either a vector crosses longitude 0 or a vector crosses
longitude -180/180. The same applies for raster maps. The general
principle here is that you must decide how to represent the earth in
latlon, and the two cases -180,180 and 0,360 should cover all
possibilities.
>
> > This would become a problem with e.g. world boundaries wrapped to
-180,180 instead of 0,360. In this case the world boundaries would need to
be transformed with x shift = 360, the original and transformed world
boundaries patched and cleaned (preferably also clipped to 0,360) and then
used as input for v.select together with the Lambert map.
>
> Unless I missed something this will work if the vector boundary does not
cross the 0/360 boundary, but will it work in a general case?
This example is devised for a vector crossing the 0,360 boundary.
Otherwise transformation without patching and cleaning would be
sufficient.
#1527: vector projection over wrapping boundary is split
-----------------------+----------------------------------------------------
Reporter: pertusus | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 6.4.2
Component: Vector | Version: svn-releasebranch64
Resolution: | Keywords:
Platform: Linux | Cpu: x86-64
-----------------------+----------------------------------------------------
Comment(by hamish):
if the raster library can handle this situation automatically and cleanly
for both the 0-360 and -180-180 cases, then why can't the vector lib be
taught do the same? forcing the user to manually do a task which should be
figured out automagically by the software seems less than ideal.
is the real problem not actually the literal numbers stored in the data
file, but rather a failure of the vector modules and display modules to
wrap LL when needed, or the vector reading library functions for not
automatically cleansing to within range?
(aka get rid of any unqualified Pythagorean shortcuts in the code, and
have libraries swap over geodesic calculations if a LL location is
detected)