polar reprojection bug revisited #2765

mapserver-bot · 2012-04-03T22:56:32Z

Reporter: gaffigan
Date: 2008/09/15 - 01:40
Trac URL: http://trac.osgeo.org/mapserver/ticket/2765
The Mapserver C library has an unresolved bug when maps requested in polar projections contain data in latlong projection. Other map/data projection combinations can also be a problem. From #723 (2004), "if the projection is valid along the borders, but contains the north pole, the source bounding box will not be correct." The result is a region of missing source data north of the incorrect northernmost northing from the source bounding box (see [attachment:ticket:2764:orig.png attached]).

Related source code:

[source:/trunk/mapserver/mapproject.c@7428#L203 mapproject.c:msProjectRect()]
[source:/trunk/mapserver/mapresample.c@7765#L1309 mapresample.c:msResampleGDALToMap()]

Previous work on this bug:

One suggestion was to initialize failure to 1 in msProjectRect(), to use a sample of points within the rectangle, instead of along the border, to compute bounds (see [http://lists.osgeo.org/pipermail/mapserver-users/2007-December/026786.html here]). In msTransformMapToSource(), setting argument bUseGrid to TRUE would be the equivalent solution.
- This requires a lot of extra operations, which are unnecessary unless the specific combination of requested projection, rectangle, and source layer projection is observed. Using a layer processing parameter would likely still result in extra operations when the request does not contain the pole.
Another suggestion was a patch to msProjectRect(), submitted in msProjectRect returns incorrect bounds resulting in features failing to show up near North or South Poles. #1999 (2006). It "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."
- The patch has not been put in place as of 5.2.0. It would need to be extended to msResampleGDALToMap(), and it does not handle cases where the source layer is not latlong.

New patch (see [attachment:ticket:2764:mapserver_polehole.patch attached]):

Affects msProjectRect() and msResampleGDALToMap.
Works by reprojecting the requested rectangle to get the source bounding box, as usual, and then repeating for an interior rectangle. If the source bounding box corresponding to the interior rectangle has a more extreme northing extent than that from the original requested rectangle, then extend source bounding box to the southernmost/northernmost extreme and the entire easting extent.
See attached [attachment:ticket:2764:ms_poletest.tar.bz2 example] and images [attachment:ticket:2764:orig.png before] and [attachment:ticket:2764:patched.png after] the patch is applied.

mapserver-bot · 2012-04-03T22:56:33Z

Author: gaffigan
Date: 2008/09/15 - 01:49
Note, I referred to attachments in the above description assuming ticket number would be 2764. Can this be changed to 2765?

mapserver-bot · 2012-04-03T22:56:34Z

Author: gosselinandre
Date: 2008/09/28 - 21:12
Replying to [ticket:2765 gaffigan]:

The Mapserver C library has an unresolved bug when maps requested in polar projections contain data in latlong projection. Other map/data projection combinations can also be a problem. From #723 (2004), "if the projection is valid along the borders, but contains the north pole, the source bounding box will not be correct." The result is a region of missing source data north of the incorrect northernmost northing from the source bounding box (see [attachment:ticket:2764:orig.png attached]).

Related source code:

[source:/trunk/mapserver/mapproject.c@7428#L203 mapproject.c:msProjectRect()]

[source:/trunk/mapserver/mapresample.c@7765#L1309 mapresample.c:msResampleGDALToMap()]

Previous work on this bug:

One suggestion was to initialize failure to 1 in msProjectRect(), to use a sample of points within the rectangle, instead of along the border, to compute bounds (see [http://lists.osgeo.org/pipermail/mapserver-users/2007-December/026786.html here]). In msTransformMapToSource(), setting argument bUseGrid to TRUE would be the equivalent solution.

This requires a lot of extra operations, which are unnecessary unless the specific combination of requested projection, rectangle, and source layer projection is observed. Using a layer processing parameter would likely still result in extra operations when the request does not contain the pole.

Another suggestion was a patch to msProjectRect(), submitted in msProjectRect returns incorrect bounds resulting in features failing to show up near North or South Poles. #1999 (2006). It "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."

The patch has not been put in place as of 5.2.0. It would need to be extended to msResampleGDALToMap(), and it does not handle cases where the source layer is not latlong.

New patch (see [attachment:ticket:2764:mapserver_polehole.patch attached]):

Affects msProjectRect() and msResampleGDALToMap.

Works by reprojecting the requested rectangle to get the source bounding box, as usual, and then repeating for an interior rectangle. If the source bounding box corresponding to the interior rectangle has a more extreme northing extent than that from the original requested rectangle, then extend source bounding box to the southernmost/northernmost extreme and the entire easting extent.

See attached [attachment:ticket:2764:ms_poletest.tar.bz2 example] and images [attachment:ticket:2764:orig.png before] and [attachment:ticket:2764:patched.png after] the patch is applied.

I unfortunately missed your ticket when I submitted mine (#2779),which appears closely related.
I tried the suggested patch, and it solved half of the problem. Now the shapes coming from the shapefiles are properly displayed. Unfortunately, the raster data are still incorrectly masked out. I have attached file bug55-patched.png to show that.

Andre Gosselin

mapserver-bot · 2012-04-03T22:56:34Z

Author: gaffigan
Date: 2008/09/29 - 04:41
Andre,

The patch file in this ticket, mapserver_polehole.patch, contains patches to two files, mapproject.c and mapresample.c, concatenated together. I notice when I [attachment:ticket:2765:mapserver_polehole.patch view this patch through trac] it only displays the first patch, to mapproject.c, which affects vector layers. The second patch, to mapresample.c, affects rasters. Check that the patch you downloaded has the following lines

--- mapresample.c.orig  2008-09-10 18:11:03.000000000 -0800
+++ mapresample.c       2008-09-13 02:04:27.000000000 -0800

To get the whole patch, you need to scroll to the bottom of the above page, to the section "Download in other formats", and click [http://trac.osgeo.org/mapserver/attachment/ticket/2765/mapserver_polehole.patch?format=raw Original Format].

When I run your example, bug55.map, on a system which has been patched, I get an image without the missing data region from the raster layer, "global". I hope you don't mind that I updated your attached image, [attachment:ticket:2765:bug55-patched.png bug55-patched.png], with my result.

Please let me know here and/or via email (gaffigan AT sfos.uaf.edu) if this image is still incorrect or if you cannot reproduce it after making sure you've patched both mapproject.c and mapresample.c. Thanks.

Steve Gaffigan

mapserver-bot · 2012-04-03T22:56:35Z

Author: gosselinandre
Date: 2008/09/29 - 06:33
Replying to [comment:3 gaffigan]:

Andre,

The patch file in this ticket, mapserver_polehole.patch, contains patches to two files, mapproject.c and mapresample.c, concatenated together. I notice when I [attachment:ticket:2765:mapserver_polehole.patch view this patch through trac] it only displays the first patch, to mapproject.c, which affects vector layers. The second patch, to mapresample.c, affects rasters. Check that the patch you downloaded has the following lines
--- mapresample.c.orig  2008-09-10 18:11:03.000000000 -0800
+++ mapresample.c       2008-09-13 02:04:27.000000000 -0800
To get the whole patch, you need to scroll to the bottom of the above page, to the section "Download in other formats", and click [http://trac.osgeo.org/mapserver/attachment/ticket/2765/mapserver_polehole.patch?format=raw Original Format].

When I run your example, bug55.map, on a system which has been patched, I get an image without the missing data region from the raster layer, "global". I hope you don't mind that I updated your attached image, [attachment:ticket:2765:bug55-patched.png bug55-patched.png], with my result.

Please let me know here and/or via email (gaffigan AT sfos.uaf.edu) if this image is still incorrect or if you cannot reproduce it after making sure you've patched both mapproject.c and mapresample.c. Thanks.

Steve Gaffigan
Steve,

I had indeed been confused by the way Trac reported the patch, and as you guessed, got only the first part dealing with file "mapproject.c". I did as you suggested and downloaded the whole patch. Now everything is OK. This is really great. Thanks so much for you help.

Not being an expert in MapServer, I do not know if this patch could create problems elsewhere in other contexts. Do you think it will need more testing before it finds its way in a new release ?

Regards,

Andre

mapserver-bot · 2012-04-03T22:56:36Z

Author: gosselinandre
Date: 2008/11/24 - 00:21
Replying to [ticket:2765 gaffigan]:
After lot of testing, I think I found a small bug in your 'mapserver_polehole.patch'.
Inside function 'msResampleGDALToMap()' the following code modifies array 'adfDstGeoTransform' :

            memcpy( &sOrigSrcExtent, &sSrcExtent, sizeof(sSrcExtent) );
            adfDstGeoTransform[0] = adfDstGeoTransform[0] + adfDstGeoTransform[1];
            adfDstGeoTransform[3] = adfDstGeoTransform[3] + adfDstGeoTransform[5];
            bSuccess =
                msTransformMapToSource( nDstXSize-2, nDstYSize-2, adfDstGeoTransform,
                                        &(map->projection),
                                        nSrcXSize, nSrcYSize,adfInvSrcGeoTransform,
                                        &(layer->projection),
                                        &sSrcExtent, FALSE );

The array is used later in the function :

/* -------------------------------------------------------------------- */
/*      Setup transformations between our source image, and the         */
/*      target map image.                                               */
/* -------------------------------------------------------------------- */
    pTCBData = msInitProjTransformer( &(layer->projection),
                                      adfSrcGeoTransform,
                                      &(map->projection),
                                      adfDstGeoTransform );

I think the call to 'msInitProjTransformer()' assumes the contents of array 'adfDstGeoTransform' to be as set at function entry, which is not the case anymore. The symptom is seen as missing lines/columns of pixels at right/bottom of resulting raster.

I was able to cure the problem by simply reinitializing array 'adfDstGeoTransform' as follows :

            memcpy( &sOrigSrcExtent, &sSrcExtent, sizeof(sSrcExtent) );
            adfDstGeoTransform[0] = adfDstGeoTransform[0] + adfDstGeoTransform[1];
            adfDstGeoTransform[3] = adfDstGeoTransform[3] + adfDstGeoTransform[5];
            bSuccess =
                msTransformMapToSource( nDstXSize-2, nDstYSize-2, adfDstGeoTransform,
                                        &(map->projection),
                                        nSrcXSize, nSrcYSize,adfInvSrcGeoTransform,
                                        &(layer->projection),
                                        &sSrcExtent, FALSE );
            /* Reset this array to its original value! */
            memcpy( adfDstGeoTransform, map->gt.geotransform, sizeof(double)*6 );

Except for that, your patch works perfectly for me.

Regards,

Andre Gosselin
Maurice Lamontagne Institute
Andre.Gosselin@dfo-mpo.gc.ca

mapserver-bot · 2012-04-03T22:56:37Z

Author: gaffigan
Date: 2008/11/25 - 09:52
Andre,

Thank you for your contribution. I've updated the attached mapserver_polehole.patch to include the reinitialization of adfDstGeoTransform.

Another possible issue with the patch is the underlying assumption: ''"If the source bounding box corresponding to the interior rectangle has a more extreme northing extent than that from the original requested rectangle, then extend source bounding box to the southernmost/northernmost extreme and the entire easting extent"''. It seems valid for the cases I work with: Combinations of map and layers in geographic, Albers, polar stereographic, and equidistant cylindrical. But there are a lot of other projections, and also rotated grids, for which the assumption needs to be evaluated.

If and when Mapserver developers get the interest/funding/time to cover this bug, a different solution may likely be preferred. Until then, I will update this ticket with any fixes or limitations I discover, and I appreciate that you're doing the same.

Steve

mapserver-bot · 2012-04-05T09:27:15Z

attachment http://trac.osgeo.org/mapserver/attachment/ticket/2765/ms_poletest.tar.bz2 :

   example map file and data

mapserver-bot · 2012-04-05T09:27:16Z

attachment http://trac.osgeo.org/mapserver/attachment/ticket/2765/orig.png :

mapserver-bot · 2012-04-05T09:27:17Z

attachment http://trac.osgeo.org/mapserver/attachment/ticket/2765/patched.png :

mapserver-bot · 2012-04-05T09:27:34Z

attachment http://trac.osgeo.org/mapserver/attachment/ticket/2765/bug55-patched.png :

   bug55-map output after applying patch (see ticket #2779)

mapserver-bot · 2012-04-05T09:27:35Z

attachment http://trac.osgeo.org/mapserver/attachment/ticket/2765/bug55-patched.2.png :

   Could not replace existing attachment.  This is the bug55.map output from a system with patches to mapproject.c and mapresample.c.

mapserver-bot · 2012-04-05T09:28:02Z

attachment http://trac.osgeo.org/mapserver/attachment/ticket/2765/mapserver_polehole.patch :

   Added reinitialization of adfDstGeoTransform (A. Gosselin)

nmtoken · 2014-04-02T13:32:42Z

I seem to be having a similar problem trying to show a tiff in native epsg:3413 projection in MapServer 6.4.1.

A request like below fails to draw any image:

http://localhost/cgi-bin/NSIDCnative/ows?language=eng&&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=-3814740.62789801321923733,-4373560.12434186320751905,1928396.42917690053582191,899341.91168433986604214&CRS=EPSG:3413&WIDTH=917&HEIGHT=841&LAYERS=ArcticRegionBlueMarbleTop60-NSIDC&STYLES=&FORMAT=image/png&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE

I tested the same tiff using GeoServer and it works (see image comparison in QGIS 2.2.0).

With MapServer I get:

With GeoServer I get:

Hoping that this isn't an issue with my tiff, which I created with:

gdalwarp -t_srs "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" arctic_decimaldegree-top60.tif arctic_decimaldegree-top60-3413.tif

on this tiff: ftp://ftp.bgs.ac.uk/pubload/OneGeology/arctic_decimaldegree-top60.tif

Incidentally, I tried first to use the arctic_decimaldegree-top60.tif (which is cut to the bounds extent of epsg:3413) and reproject to epsg:3413 (from epsg:4326), but this gives a keyhole like artefact in the map for most GetMap requests like:

mapserver-bot · 2016-02-23T17:34:03Z

This is an automated comment

This issue has been closed due to lack of activity. This doesn't mean the issue is invalid, it simply got no attention within the last year. Please reopen with missing/relevant information if still valid.

Typically, issues fall in this state for one of the following reasons:

Hard, impossible or not enough information to reproduce
Missing test case
Lack of a champion with interest and/or funding to address the issue

nmtoken · 2016-02-23T18:02:36Z

As far as I know this is still an issue

Schpidi · 2016-02-23T18:07:30Z

Is there any champion for this? Or funding to help find one?

nmtoken · 2016-02-23T18:10:42Z

Maybe there will be added interest from activities coming from http://www.opengeospatial.org/pressroom/pressreleases/2347

Schpidi · 2016-02-23T18:17:42Z

Interesting, thanks for keeping us posted.

mapserver-bot mentioned this issue Apr 3, 2012

Shapes and rasters incorrectly displayed on a map with stereographic polar (EPSG:32661) projection #2779

Closed

ghost assigned sdlime Apr 5, 2012

mapserver-bot added the No Recent Activity label Feb 23, 2016

mapserver-bot closed this as completed Feb 23, 2016

Schpidi reopened this Feb 23, 2016

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

polar reprojection bug revisited #2765

polar reprojection bug revisited #2765

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

nmtoken commented Apr 2, 2014

mapserver-bot commented Feb 23, 2016

nmtoken commented Feb 23, 2016

Schpidi commented Feb 23, 2016

nmtoken commented Feb 23, 2016

Schpidi commented Feb 23, 2016

polar reprojection bug revisited #2765

polar reprojection bug revisited #2765

Comments

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 3, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

mapserver-bot commented Apr 5, 2012

nmtoken commented Apr 2, 2014

mapserver-bot commented Feb 23, 2016

This is an automated comment

nmtoken commented Feb 23, 2016

Schpidi commented Feb 23, 2016

nmtoken commented Feb 23, 2016

Schpidi commented Feb 23, 2016