MapServer icon indicating copy to clipboard operation
MapServer copied to clipboard

Computed CRS EXTENT for EPSG 32661 UPS North - ignore axis orientation - X/Y inversion issue

Open epifanio opened this issue 9 months ago • 9 comments

Hi,

I am trying to upgrade from an older WMS service running mapserver 7.6.4 to a newer machine running mapserver 8.5 (full version below [1] [2])

  • Old service: mapserv -v [1] MapServer version 7.6.4 OUTPUT=PNG OUTPUT=JPEG OUTPUT=KML SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=SVG_SYMBOLS SUPPORTS=RSVG SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER SUPPORTS=SOS_SERVER SUPPORTS=FASTCGI SUPPORTS=THREADS SUPPORTS=GEOS SUPPORTS=POINT_Z_M SUPPORTS=PBF INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

  • New service: mapserv -v [2] MapServer version 8.5-dev PROJ version 9.4 GDAL version 3.10 OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=SVG_SYMBOLS SUPPORTS=RSVG SUPPORTS=ICONV SUPPORTS=XMP SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER SUPPORTS=OGCAPI_SERVER SUPPORTS=THREADS SUPPORTS=GEOS SUPPORTS=PBF INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE INPUT=FLATGEOBUF

The dataset used to test the two services are exactly the same and so are the mapfiles (except for path to fonts, symbols, etc)

I reduced the mapfile used on the NEW service, which is equivalent to the OLD service, to a single layer (GTiff in EPSG:32632), and looks like:

MAP
  EXTENT 5.38027653341219 60.0865467027648 7.31486520966422 61.536882894348
  # extent      5.26749797619507 60.2993491041863 7.36602770443023 61.3239944400882
  IMAGETYPE "png"
  NAME "NBS"
  SIZE 800 600
  STATUS ON
  UNITS METERS

  OUTPUTFORMAT
    NAME "png"
    MIMETYPE "image/png"
    DRIVER "AGG/PNG"
    EXTENSION "png"
    IMAGEMODE RGB
    TRANSPARENT FALSE
  END # OUTPUTFORMAT

  PROJECTION
    "init=EPSG:4326"
  END # PROJECTION
  LEGEND
    KEYSIZE 20 10
    KEYSPACING 5 5
    LABEL
      SIZE 10
      OFFSET 0 0
      SHADOWSIZE 1 1
    END # LABEL
    STATUS OFF
  END # LEGEND

  QUERYMAP
    STATUS OFF
    STYLE HILITE
  END # QUERYMAP

  SCALEBAR
    INTERVALS 4
    LABEL
      SIZE 10
      OFFSET 0 0
      SHADOWSIZE 1 1
    END # LABEL
    SIZE 200 3
    STATUS OFF
    UNITS MILES
  END # SCALEBAR

  WEB
    METADATA
      "wms_srs"    "EPSG:32632 EPSG:4326 EPSG:32661"
      "wms_onlineresource"    "https://api.get_wms_prod/f521c540-486f-42f9-9871-df3a54bf172c/wms"
      "wms_enable_request"    "*"
      "wms_abstract"    "Each of the satellites in the SENTINEL-2 mission carries a single payload: the Multi-Spectral Instrument (MSI)."
      "wms_bbox_extended"    "true"
      "wms_title"    "S2A_MSIL1C_20250120T111401_N0511_R137_T32VLN_20250120T112920"
    END # METADATA
  END # WEB

  LAYER
    DATA "false_color_glacier.tiff"
    PROJECTION
       "init=epsg:32632"
    END
    EXTENT     299995 6690245 409795 6800045
    METADATA
      "wms_srs"    "EPSG:32632 EPSG:4326 EPSG:32661"
      "wms_abstract"    "RGB Composite of bands B11, B8, B4"
      "wms_layer_group"    "/composite"
      "wms_extent"    "299995 6690245 409795 6800045"
      "wms_group_title"    "Composite"
      "wms_title"    "False Color Glacier Composite"
    END # METADATA
    NAME "false_color_glacier"

    STATUS ON
    TILEITEM "location"
    TYPE RASTER
    UNITS METERS
  END # LAYER

END # MAP

I can get the NEW service working with all the projection except for the EPSG:32661 (UPS North, which has inverted axis) when I load such layer in any client it returns blank tiles for the area where it is supposed to be - and even in QGIS, loading the layer results in a map that gets distorsion and disappears as soon as I zonn/pan on the canvas.

Trying to debug the issue, and exploring the getcapabilities of the two service, I noticed that the CRS BBOX Extent is OFF on the new service, X/Y are swapped for the projection I need to use: EPSG 32661 (which indeed as inverted axis but seems that my mapserver fails to detect that)

Of course the OLD service has no problem in let any client to render request using EPSG:32661 - see the details inb the relative getCapabilities:

<BoundingBox CRS="EPSG:32661" minx="2.2983e+06" miny="-1.35614e+06" maxx="2.43211e+06" maxy="-1.22239e+06"/>
<BoundingBox CRS="EPSG:32661" minx="-1.35552e+06" miny="2.29873e+06" maxx="-1.22275e+06" maxy="2.43151e+06"/>

Can You plerase help understanding whatr I am missing on the new service to get this to work? AS tesrting environemt on the NEW service I am using mapserver docker image from camptocamp - I built it from mapserver GIT master branch.

epifanio avatar Apr 04 '25 19:04 epifanio

I tried some more debug, I puilled an older mapserver docker image: camptocamp/mapserver:7.6-gdal3.8 and rewrote the mapfile as following:

MAP
  EXTENT 5.38027653341219 60.0865467027648 7.31486520966422 61.536882894348
  # extent      5.26749797619507 60.2993491041863 7.36602770443023 61.3239944400882
  IMAGETYPE "png"
  NAME "NBS"
  # SIZE 800 600
  SIZE    650 600
  MAXSIZE    4096
  STATUS ON
  UNITS METERS
  CONFIG "PROJ_LIB" "/usr/local/share/proj/"
  CONFIG "PROJ_DATA" "/usr/local/share/proj/"

  OUTPUTFORMAT
    NAME "png"
    MIMETYPE "image/png"
    DRIVER "AGG/PNG"
    EXTENSION "png"
    IMAGEMODE RGB
    TRANSPARENT FALSE
  END # OUTPUTFORMAT

  PROJECTION
    "init=epsg:4326"
  END # PROJECTION
  LEGEND
    KEYSIZE 20 10
    KEYSPACING 5 5
    LABEL
      SIZE 10
      OFFSET 0 0
      SHADOWSIZE 1 1
    END # LABEL
    STATUS OFF
  END # LEGEND

  QUERYMAP
    STATUS OFF
    STYLE HILITE
  END # QUERYMAP

  SCALEBAR
    INTERVALS 4
    LABEL
      SIZE 10
      OFFSET 0 0
      SHADOWSIZE 1 1
    END # LABEL
    SIZE 200 3
    STATUS OFF
    UNITS MILES
  END # SCALEBAR

  WEB
    METADATA
      "wms_srs"	"EPSG:32632 EPSG:4326 EPSG:32661"
      "wms_onlineresource"	"https://mapserver7.wps.met.no/?map=/etc/mapserver/wms_ref.map"
      "wms_enable_request"	"*"
      "wms_abstract"	"Each of the satellites in the SENTINEL-2 mission carries a single payload: the Multi-Spectral Instrument (MSI)."
      "wms_bbox_extended"	"true"
      "wms_title"	"S2A_MSIL1C_20250120T111401_N0511_R137_T32VLN_20250120T112920"
    END # METADATA
  END # WEB




  LAYER
    DATA "false_color_glacier.tiff"
    PROJECTION
       "init=epsg:32632"
    END
    EXTENT     299995 6690245 409795 6800045
    METADATA
      "wms_srs"	"EPSG:32632 EPSG:4326 EPSG:32661"
      "wms_abstract"	"RGB Composite of bands B11, B8, B4"
      "wms_layer_group"	"/composite"
      "wms_extent"	"299995 6690245 409795 6800045"
      "wms_group_title"	"Composite"
      "wms_title"	"False Color Glacier Composite"
    END # METADATA
    NAME "false_color_glacier"

    STATUS ON
    TILEITEM "location"
    TYPE RASTER
    UNITS METERS
  END # LAYER

END # MAP

I added:

  CONFIG "PROJ_LIB" "/usr/share/proj/"
  CONFIG "PROJ_DATA" "/usr/share/proj/"

and renamed the projection string to lower case:

  • at map level:
  PROJECTION
    "init=epsg:4326"
  END # PROJECTION
  • at layer level:
    PROJECTION
       "init=epsg:32632"
    END

the getcapabilities for this version+mapfile returns the expected result.

<BoundingBox CRS="EPSG:32661" minx="2.30231e+06" miny="-1.38065e+06" maxx="2.43234e+06" maxy="-1.19785e+06"/>

can you help me undertanding what I am doing wrong with mapserver 8.x ? Am I using a misconfigured mapfile or missing bits in the mapserver configuration?

Thanks!!!

EDIT: I updated to a more recent version of ubuntu/gdal while keeping mapserver to 7.6 (same issue)

epifanio avatar Apr 05 '25 08:04 epifanio

EPSG:32661 (Arctic) is marked as swapaxis but not neu in proj4 definition since a while back (unsure exactly what version of proj it did change)

also a fix for qgis have been done for the same ting

https://github.com/qgis/QGIS/pull/60597

kungzlatan avatar Apr 08 '25 08:04 kungzlatan

I have updated the service running mapserver 7.6 to include python mapscript (which I use to generate and dispatch the mapobject in my use case)

https://mapserver7py.wps.met.no/?map=/etc/mapserver/wms_ref.map&service=WMS&request=GetCapabilities

Any thought on how to proceed? As at the moment at work they require to use 32661 and I am now forced to port my app back to mapserver7.x if the issue is not addressed upstream in Mapserver

@kungzlatan should I try a more recent build of QGIS? I think the issue is on mapserver side as the older version 7.x is rendering correctly - I have been trying also client like OpenLayers - It has the same problem.

epifanio avatar Apr 08 '25 11:04 epifanio

Confirmed to be an issue/regression introduced with mapserver 8.x

From the PROJ developers,

I can't find any information abut axis orientation

really ;-) ?

    CS[Cartesian,2],
        AXIS["northing (N)",south,
            MERIDIAN[180,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",south,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[2],
            LENGTHUNIT["metre",1]],

This is a northing, easting ordered CRS. This is a MapServer issue not a PROJ one.

epifanio avatar Apr 10 '25 07:04 epifanio

@epifanio

This change that was introduced in 8.2. https://github.com/MapServer/MapServer/pull/7096. So the issue was that in earlier versions of MapServer, EPSG:32661 was wrongly treated as having axis=enu.

Regarding the GetCapabilities, I think it is correct as it is in the example for version 8.5 that you provided. According to the WMS specification, the x in <BoundingBox> corresponds to the first axis in the CRS (i.e. y in 32661). See discussion at: https://github.com/MapServer/MapServer/issues/6978

So if we transform the minx and miny coordinates from the BBOX in you example, I get values close to the extent defined in the MAP file:

<BoundingBox CRS="EPSG:32661" minx="-1.35552e+06" miny="2.29873e+06"

minx=-1.35552e+06
miny=2.29873e+06

x=$miny # x in <BoundingBox> corresponds to y in 32661
y=$minx # y in <BoundingBox> corresponds to x in 32661
echo $x $y | gdaltransform -s_srs EPSG:32661 -t_srs EPSG:4326 # 5.08742708186411 60.3127167369556 0

PatrikSylve avatar Apr 28 '25 12:04 PatrikSylve

The GetCapabilities links do not work at the moment, but next time it would be good to write the request with all parameters, expecially including &VERSION=1.3.0 because the axis order thing was changed between versions 1.1.1 and 1.3.0 (now somebody will say that it was not changed but it was only understood in a wrong way before - it means the same from users perspective).

Maybe Mapserver is using &version=1.3.0 by default but it does not hurt to be explicit.

jratike80 avatar Apr 28 '25 15:04 jratike80

yes, I was forcing the version to be 1.3 - I create the getcapabilities from inside FastAPI via mapscript.

As suggested in https://github.com/MapServer/MapServer/issues/7273 - I rebuilt QGIS using git master - and it is now loading the WMS layer directly from the endpoint even when specifying the CRS EPSG:32661

This is great, as it seems that now both QGIS / Mapserver / GDAL / Proj are aligned in using the same axis convention when dealing with inverted axis projection.

As mentioned in https://github.com/MapServer/MapServer/issues/7273 the problem will now arise when loading the WMS into a JS client - to see if it can handle EPSG:32661 correctly.

I have some OL code to work wit proj-js and switch between EPSG:4326 and EPSG:32661 but this was using mapserver7.x - the same client code fails to render the WMS when provided by mapserver V 8.x

I will try to make sure I have the latest version of both OL and proj-js and try again - and report back here @PatrikSylve @jratike80 - thanks a lot for looking itno this :)

epifanio avatar Apr 28 '25 17:04 epifanio

this is now getting interesting, I changed the FastAPI mapserver API to use a jinja template and create a mapfile (based on mapserver 7.x ) which is then served via mapscript - this is using mapscript + mapserver V7.x

A temporary endpoint: https://map.wps.met.no/get_wms/26c3a197-c590-4e2b-bd44-ff213b56a1ce/wms

this works as expected in the "old" QGIS - I tried V3.34.4-Prizren - but not in QGIS master (I guess this was gonna be expected)

I guess my only alternative to have a stable service that will work with any client is to implement some logic to detect the client that is sending the request ..

I may try to build 2 jinja templates (mapserver V7 + V8) and instead of dispatching the service via FastAPI+mapscript I will have to redirect the response to a mapserver fastcgi service. like having 2 mapserver fcgi services running in different containers (one for V7 and an other for V8) and let the FastAPI redirect to the mapserv version based on some information from the client request :( it seems painful :(

epifanio avatar Apr 29 '25 12:04 epifanio

Hi,

I have QGIS 3.43.0, installed with OSGeo4W. I set the QGIS project CRS into EPSG:32661 and asked WMS to use the same CRS. Zoom to layer works OK but zooming in returns results that look unexpected to me and they do not fill the bbox that I zoomed in. But maybe it is normal because the uncommon asis orientation of this coordinate system https://geohubkenya.wordpress.com/2019/08/15/ups-coordinate-system/

One request captured from QGIS as an example.

https://map.wps.met.no/get_wms/26c3a197-c590-4e2b-bd44-ff213b56a1ce/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=1120010%2C-692280%2C1702890%2C90152.19999999999709&CRS=EPSG%3A32661&WIDTH=159&HEIGHT=120&LAYERS=composite&STYLES=&FORMAT=image%2Fpng&DPI=144&MAP_RESOLUTION=144&FORMAT_OPTIONS=dpi%3A144&TRANSPARENT=TRUE

jratike80 avatar Apr 29 '25 13:04 jratike80