Opened 12 years ago

Last modified 12 years ago

#2765 new defect

polar reprojection bug revisited

Reported by: gaffigan Owned by: sdlime
Priority: normal Milestone:
Component: MapServer C Library Version: unspecified
Severity: normal Keywords:
Cc: warmerdam

Description

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:

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.

Attachments (6)

ms_poletest.tar.bz2 (29.3 KB) - added by gaffigan 12 years ago.
example map file and data
orig.png (45.4 KB) - added by gaffigan 12 years ago.
patched.png (75.0 KB) - added by gaffigan 12 years ago.
bug55-patched.png (114.5 KB) - added by gosselinandre 12 years ago.
bug55-map output after applying patch (see ticket #2779)
bug55-patched.2.png (192.1 KB) - added by gaffigan 12 years ago.
Could not replace existing attachment. This is the bug55.map output from a system with patches to mapproject.c and mapresample.c.
mapserver_polehole.patch (6.7 KB) - added by gaffigan 12 years ago.
Added reinitialization of adfDstGeoTransform (A. Gosselin)

Download all attachments as: .zip

Change History (13)

Changed 12 years ago by gaffigan

Attachment: ms_poletest.tar.bz2 added

example map file and data

Changed 12 years ago by gaffigan

Attachment: orig.png added

Changed 12 years ago by gaffigan

Attachment: patched.png added

comment:1 Changed 12 years ago by gaffigan

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

comment:2 in reply to:  description Changed 12 years ago by 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:

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

Changed 12 years ago by gosselinandre

Attachment: bug55-patched.png added

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

comment:3 Changed 12 years ago by 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

Changed 12 years ago by gaffigan

Attachment: bug55-patched.2.png added

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

comment:4 in reply to:  3 Changed 12 years ago by 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

comment:5 in reply to:  description Changed 12 years ago by 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@…

Changed 12 years ago by gaffigan

Attachment: mapserver_polehole.patch added

Added reinitialization of adfDstGeoTransform (A. Gosselin)

comment:6 Changed 12 years ago by 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

comment:7 Changed 12 years ago by warmerdam

Cc: warmerdam added
Note: See TracTickets for help on using tickets.