Quebec separated from Canada (!!!)

[dmorissette]

This is not a political post... but the above WMS URL does not seem consider 
Quebec to be part of Canada any more! :)

Actually, what happens is that with an image size of 300x300, Quebec and the 
Maritimes are gone, but if you change the size to 600x300 then all the provinces 
are fine.  Interesting!?!??!?!

[dmorissette]

Image produced by this bug's URL and showing the problem

[jdoyon@…]

I just noticed this on my end as well (Daniel and I use the same data).

I'm not sure, my memory could be failing me, but I *think* this started 
happening when I updated my code from the CVS today, my previous MS3.5 was from 
a couple of weeks ago ... Maybe changes that occured in the mean time caused 
this ?

That's one possiblity, but I also noticed that they disappear only when 
reprojecting, and then, only certain projections.  If i do a very simple OGC 
call (?wmtver=1.0.0&request=map) with any size, the missing provinces are 
allways there. When I reproject using an OGC client using EPSG:4326, then 
Quebec and Maritimes are missing, BUT, if I reproject to EPSG:42101 (WSG 84 / 
LCC Canada), then the provinces show up fine ...

Hope this helps...

[fwarmerdam]

In maputil.c:msDrawLayer() (and presumably it's analogs) the following
call is done to transform the latlong region of the map into the projection
of the layer (LCC).  I think this transformation is operating very poorly
because the lat/long region is so unreasonable as input to the LCC reprojection.

  if((map->projection.numargs > 0) && (layer->projection.numargs > 0))
    msProjectRect(&map->projection, &layer->projection, &searchrect); // projec

(gdb) p map->extent
$19 = {minx = -179.40000000000001, miny = -179.40000000000001, 
  maxx = 179.40000000000001, maxy = 179.40000000000001}
(gdb) p searchrect
$20 = {minx = -2306170961.9598298, miny = -566618157.16840172, 
  maxx = 3.1311206780778273, maxy = 4642261.9173844298}
(gdb) 

Some shapes are failnng the subsequent test against searchrect, including 
Quebec, and Greenland.



On review of msProjectRect(), a good deal of effort is expended to 
compute a bounding rectangle that is likely to be OK.  A sampling
of points around the edge are used, not just the corners. 

However, it doesn't account for the fact that pj_transform() can fail
for illegal points, and so the bounds are not well computed.  All the top and
bottom edges have no useful effect (but are misconstrued as being small values
since the failure of pj_transform() is not detected).  The right side points
are projected to have large negative x values which is reasonable.  The
problem is that we don't do points through the center that can approach the
right edge in a valid region, getting a properly MAXX value. 

One possible solution is to to the current sampling all through the rect
not just along the edges but this will substantially increase the computational
cost.  Perhaps this would just kick in when a transform fails. 

[fwarmerdam]

I have fixed the bug by improving msProjectRect()'s ability in the face
of some points failing to reproject. 

msProjectPoint now returns MS_FAILURE if transform fails.
msProjectRect now deals with failure to reproject some points, and samples
through the middle of the rectangle if there are any failures.