Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Incorrectly rendered EPSG:4326 raster in Equal Earth, depending on output image size #5958

Open
dbaston opened this issue Dec 31, 2019 · 7 comments
Milestone

Comments

@dbaston
Copy link
Contributor

@dbaston dbaston commented Dec 31, 2019

When rendering a raster stored in EPSG:4326 to EPSG:8857 (Equal Earth), values are plotted in incorrect locations for certain ouput image sizes.

Mapfile:

MAP
  PROJECTION 
    "init=epsg:4326"
  END

  EXTENT -180 -90 180 90        

WEB
  METADATA
  "ows_enable_request" "*"
  "wms_srs"            "EPSG:4326 EPSG:8857" 
  END
END

LAYER
  NAME "rast"
  TYPE RASTER
  STATUS OFF 

  PROJECTION "init=epsg:4326" END

  DATA "/tmp/data.tif"

  CLASS
    EXPRESSION ([pixel] >= 0)
    STYLE
      COLORRANGE 0 0 0 0 0 255
      DATARANGE 0 60
      RANGEITEM "pixel"
    END
  END
END

END

Raster data: data.zip

Expected output:

./mapserv -nh "QUERY_STRING=map=/home/dan/dev/mapserver/test.map&service=wms&version=1.3.0&request=GetMap&format=png&width=2162&height=1055&crs=EPSG%3A8857&layers=rast&BBOX=-17243959.06221695,-8392927.59846645,17243959.06221695,8392927.59846645&transparent=TRUE" > ok.png

image

Incorrect output (image width increased from 2162 to 2163 pixels):

./mapserv -nh "QUERY_STRING=map=/home/dan/dev/mapserver/test.map&service=wms&version=1.3.0&request=GetMap&format=png&width=2163&height=1055&crs=EPSG%3A8857&layers=rast&BBOX=-17243959.06221695,-8392927.59846645,17243959.06221695,8392927.59846645&transparent=TRUE" > bad.png

image

./mapserv -v
# MapServer version 7.5-dev OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WFS_SERVER SUPPORTS=WCS_SERVER SUPPORTS=FASTCGI SUPPORTS=GEOS SUPPORTS=PBF INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE
@dbaston

This comment has been minimized.

Copy link
Contributor Author

@dbaston dbaston commented Dec 31, 2019

The two cases (width=2162 or width=2163) and up hitting a different path around here:
https://github.com/mapserver/mapserver/blob/master/mapresample.c#L976

In the 2162 case, dfDeltaX is computed to be positive, and the resampling takes place by recursing several times into subsets of the scanline. In the 2163 case dfDeltaX is computed to be negative and the resampling takes place with a single pass on the scanline.

@dbaston

This comment has been minimized.

Copy link
Contributor Author

@dbaston dbaston commented Dec 31, 2019

Here's what seems to be going on.

  • msApproxTransformer wants to transform a scanline in map pixel space (0.5 to 2162.5) to image space (0-719).
  • It first uses PROJ to transform the first, middle, and last points in pixel space (0.5, 1081.5, and 2162.5) to image space and gets (719.73, 360, 0.267)
  • It then tries to use linear interpolation from the first and last points to approximate the middle point. The intention is to avoid using PROJ to transform all of the points if linear interpolation is accurate enough.
  • In this case, the linear interpolation exactly computes the middle point but is wrong for all other points! If we modify the code so that entire scanline is transformed by PROJ, we get values that look like this: 719.73, 0.065, 0.398, ... 719.26, 219.60, 719.93, 0.267. So the delta estimated by linear interpolation from the first and last points has the wrong sign, and all intermediate points (except the middle one!) are incorrect. And this only happens when the first and last points in the scanline are extreme enough, which is why we only see it for output image widths >= 2163 pixels.

A workaround is to disable linear interpolation by initializing the approximate transformer with a negative acceptable error on line 1768.

@dbaston

This comment has been minimized.

Copy link
Contributor Author

@dbaston dbaston commented Dec 31, 2019

This particular case could be fixed by using PROJ to transform five points instead of three and then checking that, in addition to a successful interpolation of the center point, the sequence of test points is monotonic. But I'm not sure how robust of a fix that is.

@dbaston

This comment has been minimized.

Copy link
Contributor Author

@dbaston dbaston commented Jan 3, 2020

I think this was closed by mistake, @rouault ?

@rouault

This comment has been minimized.

Copy link
Contributor

@rouault rouault commented Jan 3, 2020

I think this was closed by mistake, @rouault ?

yes, sorry. I've used a wrong ticket number in the commit message of a unrelated fix. Re-opening

@rouault rouault reopened this Jan 3, 2020
@rouault

This comment has been minimized.

Copy link
Contributor

@rouault rouault commented Jan 3, 2020

@dbaston Great analysis of this issue by the way. Perhaps a potential solution would be not to use exactly the middle point, but a point at 45% of the segment . While not being fully bullet proof, it could probably be good enough.

dbaston added a commit to dbaston/mapserver that referenced this issue Jan 7, 2020
This corrects for an edge case where linear interpolation is successful
for the middle point of a scanline but is incorrect for all other
points.

Fixes mapserver#5958
dbaston added a commit to dbaston/mapserver that referenced this issue Jan 7, 2020
This corrects for an edge case where linear interpolation is successful
for the middle point of a scanline but is incorrect for all other
points.

Fixes mapserver#5958
@jmckenna jmckenna closed this in 5f89b67 Mar 22, 2020
@dbaston

This comment has been minimized.

Copy link
Contributor Author

@dbaston dbaston commented Mar 24, 2020

@jmckenna can you please re-open? The commit message in 5f89b67 references this issue by accident.

@sdlime sdlime reopened this Mar 24, 2020
@jmckenna jmckenna added this to the 8.0 Release milestone Mar 24, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked pull requests

Successfully merging a pull request may close this issue.

None yet
4 participants
You can’t perform that action at this time.