[MapServer-users] Contour and wind barbs problems with later MS, GDAL, Proj versions

[David Miller]

 Good day Mapserver Users!

Not too long ago, we upgraded our development area to use Mapserver 7.6.4
which linkage to GDAL 3.4.0 and Proj 7.2.  Unfortunately, the contour
images created via MS over areas which span the dateline are now wrapping
back from one end to the other as shown in the first image at
https://pasteboard.co/rHXaEwQEOT13.png.

Our production area is using older versions (MS 7.2.2, GDAL 2.4, and Proj
4.8) does not have this problem as indicated at
https://pasteboard.co/aJL1HE6lIZ1u.png.

This issue occurs with at least two of our datasets where the rasters which
span the dateline and where we're using the MS CONNECTIONTYPE CONTOUR (and
as I note later, this affects wind vector point creation too,
CONNECTIONTYPE UVRASTER).  The longitudes in the images above span from
120E to 1W.  In order for MS to produce this image (originally from GRIB
files) and also put it the Mercator projection we wish to display it in
(mostly because such rasters from our other datasets are in Mercator
natively), we warp the grid from lat/lon using gdalwarp:

gdalwarp -q -s_srs '+proj=latlong +lon_wrap=180 +datum=WGS84 +no_defs
+unit=degree +a=6371229 +b=6371229' -t_srs "+proj=merc +a=6378137
+b=6378137 +lat_0=0.0 +lat_ts=0.0 +lon_0=180.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m" -ts 1912 0 -ot Int16 -overwrite -co COMPRESS=DEFLATE -co
BLOCKXSIZE=128 -co BLOCKYSIZE=128 -co ZLEVEL=9 -co TILED=YES -co PREDICTOR=2


As indicated, the -s_srs proj string is using the +lon_wrap=180 to indicate
the longitudinal wrap and the -t_srs proj string is indicating lon_0 as
180, also indicating the need to produce a dateline crossing image.

In the MS epsg file, we created this projection for the contour image
returned from MS:

<920919> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=-180.0
+x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs +over <>

The map projection overall is

<920916> +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +k=1
+units=m +no_defs +over <>

This is used for the data raster overlay in PNG format.

In the mapfile, these would appear as such (I included a bit more of the
LAYER definition):

In the MAP section:

PROJECTION
       "init=epsg:920916"
  END

In the LAYER section:

LAYER
    NAME gcntr
    TYPE LINE
    STATUS DEFAULT

    CONNECTIONTYPE CONTOUR
....

    PROJECTION
   "init=epsg:920919"
    END

    PROCESSING "BANDS=1"
    PROCESSING "CONTOUR_ITEM=cvalue"
    PROCESSING "CONTOUR_INTERVAL=%intval%"
    PROCESSING "CONTOUR_LEVELS=%levels%"
    GEOMTRANSFORM (smoothsia(generalize([shape], 0.25*[data_cellsize])))
    LABELITEM "cvalue"
..... (EXPRESSION statements follow)

  The +lon_0=-180.0 for the contour layer projection properly positions the
contours on top of the image overlay using the MAP's projection like so
using Openlayers 4.6.5 on our production server as shown in
https://pasteboard.co/DTSlmj4cBt0V.png.

Same mapfile with the same epsg file on our development platform with the
updated MS, GDAL, and Proj yields the following display when rendered via
Openlayers, https://pasteboard.co/YzryouudsshD.png.

The contours are not only wrapping but are cut off west of 180.
Interestingly, the color dataset PNG image displays properly.  If we pan
the map in our viewer where the center is east of the dateline, the
horizontal wraparound line disappears.  Needless to say, we're very
reluctant to update our production MS, GDAL and Proj with this issue
present.  It also affects wind vector images as well BTW.

Just for more information, we're slicing out global-spanning GRIB datasets
to the window I indicated, i.e. on a 0-360 dataset, it goes from 120 to 359
(or 120E to 1W) and 80N to 30S (  we had run into problems initially going
to the Prime Meridian).  We process the subset area in GDAL converting
units and also saving to GeoTIFF as those have seemed to be easier working
wrt geospatial considerations.  And finally, that GeoTiff is warped from
the GRIB data's native projection to one we wish to display it in, in this
case Mercator as a couple of other different dateline spanning datasets
contain the data in that projection.

We have tried different combinations of not using the +lon_wrap as well as
+lon_0=-180 and also +over.  However, we always end up with the contour
lines wrapping around and the lines cut off to the west of the dateline.
We have also experimented with using different combinations of MS w/GDAL
and Proj. old and new:

MS 7.2.2, GDAL 3.4, Proj 7.2 no contour wrap
MS 7.2.2, GDAL 3.4, Proj 4.8 no contour wrap
MS 7.6.4, GDAL 3.4, Proj 4.8 contour wrap

And of course the original configurations

MS 7.2.2, GDAL 2.2, Proj 4.8 no contour wrap (production server)
MS 7.6.4, GDAL 3.4, Proj 7.2 contour wrap (development server)

I searched release notes for wrapping fixes and only a few were listed in
7.4.4 - #5975, 5965, 5960.   Perhaps one of these accidentally affected the
MS CONTOUR and UVRASTER functions for dateline-spanning rasters?

I'll create a bug report if needed but thought I'd ask this group first.

Thanks for any help.

David Miller
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/mapserver-users/attachments/20220830/58584895/attachment.htm>