polar reprojection bug revisited

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

Related source code:

• mapproject.c:msProjectRect()​
• 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 ​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 #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 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 example and images before and after the patch is applied.

[gaffigan]

example map file and data

[gaffigan]

Note, I referred to attachments in the above description assuming ticket number would be 2764. Can this be changed to 2765?

[gosselinandre]

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

Related source code:

• mapproject.c:msProjectRect()​
• 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 ​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 #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 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 example and images before and 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

[gosselinandre]

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

[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 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 ​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, 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

[gaffigan]

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

[gosselinandre]

Replying to 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 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 ​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, 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

[gosselinandre]

Replying to 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@…

[gaffigan]

Added reinitialization of adfDstGeoTransform (A. Gosselin)

[gaffigan]

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