vector projection over wrapping boundary is split

[pertusus]

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.

[mmetz]

patch to wrap ll coords from 0 to 360 instead of from -180 to 180

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

Markus M

[pertusus]

There are those compilation warnings with -Wall.

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.

[neteler]

Perhaps these reports should remain open until the backporting is done?

[mmetz]

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.

Markus M

[pertusus]

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

I attach a tarball similar but for 6.4.

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

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

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

Markus M

[pertusus]

The projection indeed seems to do the same than with the previous patch. The other issues are still there, though...

[mmetz]

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

Markus M

[pertusus]

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

ok.

[mmetz]

Replying to pertusus:

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

Markus M

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

?, Hamish (170E)