msProjectRect returns incorrect bounds resulting in features failing to show up near North or South Poles.

[sroberts@…]

The subroutine msProjectRect in mapproject.c returns incorrect bounds when the
output bounds are in latitude/longitude, the input projection is something other
than latlon and either the North or South Pole lies inside the input bounds. One
symptom of this error is when using a standard projected map where the pole is
in view, features specified by latitudes/longitudes that are near the pole will
not show up.

Below is a url to an example mapserver site where I am trying to plot ship
tracks and coring sites that cross the North Pole on a standard Polar
Stereographic map:

http://mapserver.eol.ucar.edu/cgi-bin/bug/mapserv-4.10.0?mode=browse&map=/export/web/mapserver/data/bug/test.map&imgxy=300+300&imgext=-2000000+-2000000+2000000+2000000&zoomsize=2

I'm using mapserver 4.10.0 in this example. As you can see not much shows up
with the initial map extent which is centered on the North Pole. But when you
recenter the map away from the pole more of the ship tracks and core sites will
show up.

The core of the problem is based on the false assumption that the lat/long
equivalent of the edge of the map bounding box can be used to represent a
rectangle bounding box to filter lat/long data that is outside of your map view.
In the above example url with the map centered on the North Pole the edges of
the bounding box never gets above 71.5 degrees so we end up with a circular
filter "box" of 180W,180E,64N,71.5N where everything outside of this "box" is
not plotted. Thus everything north of 71.5N is not plotted! This is obviously
not what we want. There is an excellent article on the web which discusses these
issues at:

  http://chukchi.colorado.edu/PAPER/
  
msProjectRect includes logic to "sample" through the destination rectangle if
there is a failure to reproject any of the edge points. In this example and many
others, reproject failures never occur so this does not solve the problem. And
doing a course sampling through the interior of the rectangle is not an optimal
solution either. There is a high probability of missing the pole. In my example
url above I have a coring site right at the North Pole and when I modify the
code to force "sampling" most of the data will now show up but the North Pole
core site still fails to plot. A proposed simple solution would be to add some
code to detect if the poles are inside the map bounding box and adjust the
extents accordingly. Adding the following code to the subroutine msProjectRect
in mapproject.c just before the line "if( failure > 0 )" would do the trick:

   /* Check if North or South Pole is in search rectangle */
   if(!pj_is_latlong(in->proj) && pj_is_latlong(out->proj) ) {
       /* Check if North Pole is in search rectangle */
       prj_point.x = 0;
       prj_point.y = 90;
       if( msProjectPoint( out, in, &prj_point ) == MS_FAILURE ) {
           msDebug( "msProjectRect(): Failed to project North Pole.\n" );
           pole_failure++;
       } else {;
           if(msPointInRect(&prj_point,rect)) {
              prj_rect.maxy = 90;
           }
       }

       /* Check if South Pole is in search rectangle */
       prj_point.x = 0;
       prj_point.y = -90;
       if( msProjectPoint( out, in, &prj_point ) == MS_FAILURE ) {
           msDebug( "msProjectRect(): Failed to project South Pole.\n" );
           pole_failure++;
       } else {;
           if(msPointInRect(&prj_point,rect)) {
              prj_rect.miny = -90;
           }
       }
       if(pole_failure > 1) {
           msDebug( "msProjectRect(): Failed to project both North and South
Pole.\n" );
           failure++;
       }
   }


And you will need to add:

   int     pole_failure=0;

at the beginning of the subroutine. The above code works by projecting the north
and south poles from lat/longs to the map projection and then using the routine
"msPointInRect" to check if these points are inside the map bounding box. With
the above addition the above example will now generate the bounds
180W,180E,64N,90N and all the data should now show up.

However, I've noticed that msProjectRect still has issues when the lat/long
bounding box crosses the dateline. Regardless of how small the bounding box is,
if it crosses the dateline it will produce a longitudinal extent of -180 to 180.
 This will potentially generate a large number of "False Positives" and possibly
cause other more serious problems.

-Steve

[sdlime]

Frank's the projection guru so cc'ing...

[sroberts@…]

Contains a map file and sample shape files that demonstrate this bug.