Vector Reprojection problem with whole world extent (Alaska)

[dmorissette]

With the above URL using the GMap data and trying to display the whole world 
extent, some points from Alaska are reprojected at the other end of the world 
and this creates lines across the top of the map.

We need to find a way to at least avoid this (clipping?), and if possible 
properly draw the portion of Alaska that should go at the top-right of the map.

[dmorissette]

The map produced by this bug's URL and showing the bug

[fwarmerdam]

I have setup the data on my system, and the problem can be seen in this
request:

http://gdal.velocet.ca/cgi-bin/mswms_gmap?wmtver=1.0.0&request=map&BBOX=-180,-90,180,90&WIDTH=600&HEIGHT=300&SRS=EPSG:4326&LAYERS=land_fn

PROJ.4 is automatically wrapping longitude values in the range -180 to 180 so
a line (reprojected from LCC to lat/long) that happens to cross the dateline
wraps from -180 to +180 causing the horizontal lines stretched across the world.


One approach to resolving this would be to ensure that lat/long values are
not wrapped.  However, this would require scary internal changes to PROJ.4, and
it isn't clear to me that this is a reasonable thing to do, partly because
it puts a magical importance on the dateline instead of incorporating an
understanding that lat/long is special and "wraps" unlike projections. 

While investigating this, I noticed something else. MapServer won't do anything
meaningful with extents outside of -180 to 180.  For instance the following
url going from -110 to +250 (intended to be another world view with a different
cutline) doesn't draw anything in the area of +180 to +250. 

http://gdal.velocet.ca/cgi-bin/mswms_gmap?wmtver=1.0.0&request=map&BBOX=-110,-90,250,90&WIDTH=600&HEIGHT=300&SRS=EPSG:4326&LAYERS=land_fn

[fwarmerdam]

Discovered that the wrapped view can be correctly generated with the following
URL, using a BBOX MAXX value slightly smaller than the MINX value. 

http://gdal.velocet.ca/cgi-bin/mswms_gmap?wmtver=1.0.0&request=map&BBOX=-110,-90,-110.000000001,90&WIDTH=600&HEIGHT=300&SRS=EPSG:4326&LAYERS=land_fn

[fwarmerdam]

Daniel / Steve,

I have been assigned reprojection bug:

  http://mapserver.gis.umn.edu/bugs/show_bug.cgi?id=15

I have hacked together a solution in mapproject.c by assuming that
lines should take the short route around the earth and so to force
longitudes to go outside the range -180 to 180 if it will keep a line
"tight".  The result from our test case is here:

http://gdal.velocet.ca/cgi-bin/mswms_gmap?wmtver=1.0.0&request=map&BBOX=-180,-90,180,90&WIDTH=600&HEIGHT=300&SRS=EPSG:4326&LAYERS=land_fn

Basically I modify msProjectLine() to do a simple test:

  if( be_careful )
  {
      for(i=0; i<line->numpoints; i++)
      {
          double        dist;

          msProjectPoint(in, out, &(line->point[i]));
          if( i > 0 )
          {
              dist = line->point[i].x - line->point[i-1].x;
              if( dist > 180 )
              {
                  line->point[i].x -= 360.0;
              }
              else if( dist < -180 )
              {
                  line->point[i].x += 360.0;
              }
          }
      }
  }

Of course this doesn't do the right thing if we were trying to make a
really wide box (eg from -175 to +175). 

I can apply further analytical improvements to essentially always do the
right thing by reprojecting a point in the middle of each line segment
and testing to see if we are closer to the middle of the line segment
in lat/long space if we go around the earth the other way.  However,
this test would be fairly expensive, at least doubling the amount of
reprojection to be done. 

Furthermore, a more complete solution is to use such rules to recognise
line segments that should cross the dateline and actually splitting the
line into separate segments approaching the dateline from each side. 
This would require rearranging the shapeObj during splitting. Something
similar would need to be done to split polygons along the dateline. 

The problem is that all these solutions are increasingly expensive
in CPU time, with substantial computation added for all reprojected
vectors. 

So the question is; how much do we want to do?  Do we want to control
the expensive operations on the basis of some configurable (compile time
or mapfile?) parameter? 

Best regards,

---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

[dmorissette]

Subject:             Re: Wrapping problems with Lat/long
   Date:             Fri, 06 Jul 2001 16:57:18 -0400
   From:             Frank Warmerdam <warmerdam@pobox.com>
     To:             Daniel Morissette <morissette@dmsolutions.ca>
     CC:             Steve Lime <steve.lime@dnr.state.mn.us>

Daniel Morissette wrote:
> 
> Frank,
> 
> How about this, it's a variation of your middlepoint suggestion:
> 
> At the beginning of msProjectLine(), compute the MBR of the whole
> polyline and reproject its corners, if they appear to be crossing the
> date line, then reproject the center of the MBR, and check if the MBR
> center is closer to the date line than the corners are... depending on
> the result, set a flag that will be used for all the coordinates in the
> polyline.
> 
> It should give a result similar to your middlepoint reprojection
> suggestion, except that you do the middlepoint test only once so it's
> not too expensive.

Daniel, 

This would indeed work, and help for complex polylines/polygons.

> BTW, after thinking about it more, I'm surprised that your little test
> worked... what happens with line segmetns completely WEST of the date
> line but that are part of a polygon that crosses the date line?  Your
> test should catch only the line segments that do cross the date line, so
> I'm surprised that the URL you gave us did work... unless it was a lucky
> coincidence and there was only one single point in the whole polygon
> that was on the WEST side of the date line...

The trick is that you fix up each point before testing the next.  So 
basically all the points after the first get dragged to the same side
of the world as the first.  It's induction!

I am going to put off fixing the wrapping problem till tomorrow to give
us all time to think about it a bit.  However, I did bring Quebec back 
into Canada.  That was a sticky but solvable problem (at least within
MapServer). 

[fwarmerdam]

I have reviewed the suggestion related to using the MBR and the center point,
but I can't work this through in such a way that it is safe and always works. 
Instead I have come up with this proposal.


Proposal:

Modify msProjectLine() so that it "dateline wraps" objects when necessary
in order to preserve their shape when reprojecting to lat/long.  This
will be accomplished by:

1) As each vertex is reprojected, compare the X distance between that 
   vertex and the previous vertex.  If it is less than 180 then proceed to
   the next vertex without any special logic, otherwise:

2) Reproject the center point of the line segment from the last vertex to
   the current vertex into lat/long.  Does it's longitude lie between the
   longitudes of the start and end point.  If yes, return to step 1) for
   the next vertex ... everything is fine. 

3) We have determined that this line segment is suffering from 360 degree
   wrap to keep in the -180 to +180 range.  Now add or subtract 360 degrees
   as determined from the original sign of the distances.  

This is similar to the code there now (though disabled in CVS); however, 
it will ensure that big boxes will remain big, and not get dateline wrapped
because of the extra test in step 2).  However step 2 is invoked only very
rarely so this process takes little more than the normal process.  In fact, 
if we were sufficiently concerned about performance we could do a test on the
shape MBR in lat/long space, and if the width is less than 180 we know we never
need to even do test 1). 

What doesn't this resolve:

This ensures that individual lines are kept in the proper shape when 
reprojected to geographic space.  However, it does not:

 o Ensure that all rings of a polygon will get transformed to the same "side"
   of the world.  Depending on starting points of the different rings it is
   entirely possible for one ring to end up in the -180 area and another ring
   from the same polygon to end up in the +180 area.  We might possibly be
   able to achieve this though, by treating the multi-ring polygon as a whole
   and testing the first point of each additional ring against the last
   vertex of the previous ring (or any previous vertex for that matter).

 o It does not address the need to cut up lines and polygons into distinct
   chunks to preserve the correct semantics.  Really a polygon that 
   spaces the dateline in a -180 to 180 view should get split into two 
   polygons.  We haven't addressed that, though if it were to be addressed,
   it could be done as a followon and distinct step from what we are doing
   above.  In fact this sort of improvement (split polygons based on dateline
   or view window) should be done for all lat/long shapes regardless of 
   whether they are being reprojected from another projection. 

 o It does not address issues related to viewing rectangles that go outside
   the -180 to 180 longitude range.  For instance, it is entirely reasonable
   to want a 160 to 200 longitude view to see an area on the dateline clearly.
   Currently shapes in the -180 to -160 range which should be displayed in the
   180 to 200 portion of that view will not be because there is no recogition
   that they belong there. 

[fwarmerdam]

Proposed Changes applied to mapproject.c. I have tested by introducing an
extra large rectangle, and verified that it was not "wrapped" while the
siberia area was properly wrapped.  

bug fixed.