MapServer icon indicating copy to clipboard operation
MapServer copied to clipboard

Published 20 hours ago • verified publisher icon MapServer

CONNECTIONTYPE CONTOUR + TILEINDEX gives weird results

Open landryb opened this issue 8 years ago • 7 comments

On 7.0.4, i'm trying the awesome CONTOUR feature, and i've had mixed results with this LAYER definition:

LAYER
        NAME "rgealti_5m_courbes_niveau"
        MAXSCALEDENOM 20001
        INCLUDE "common/wfs_queryable.map"
#       TILEINDEX "dallage_rgealti_5m_mnt"
        TILEINDEX "/data/wxs/ign/rgealti/2016/MNT/dallage.shp"
        TILEITEM "location"
#       DATA "/data/wxs/ign/rgealti/2016/MNT/all.vrt"
        TYPE LINE
        CONNECTIONTYPE CONTOUR
        PROCESSING "CONTOUR_ITEM=elevation"
        PROCESSING "CONTOUR_INTERVAL=0,5000:10"
        PROCESSING "CONTOUR_INTERVAL=5000,10000:50"
        PROCESSING "CONTOUR_INTERVAL=10000,20000:100"
        GEOMTRANSFORM (smoothsia([shape], 5))
        LABELITEM "elevation"
        CLASSITEM "elevation"
....
....
END

If i use one of both TILEINDEX sources (the first one references another layer in the same mapfile, which is a vector layer with a postgis connection), the contour lines are only computed on a single source tile and not across all tiles covered by the bbox query - more or less not always the one at the center of the bbox extent but one on the side of the bbox border. See on the screenshot, the black line is the tileindex rendering. contour

Of course no problem if i use the DATA source pointing at the vrt referencing the same source tiles.

landryb avatar Feb 15 '17 16:02 landryb

Using shp2img, here's what i have in the debug log for the failing case, pointing TILEINDEX at the shapefile:

[Wed Feb 15 17:19:41 2017].461991 Entering msContourLayerWhichShapes().
[Wed Feb 15 17:19:41 2017].461997 msConnPoolRelease(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7f45fda2ca00)
[Wed Feb 15 17:19:41 2017].462002 msOGRLayerClose(__rgealti_5m_courbes_niveau_CONTOUR__).
[Wed Feb 15 17:19:41 2017].462005 msOGRFileClose(,0).
[Wed Feb 15 17:19:41 2017].462009 msConnPoolRelease(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7f45fda2ca00)
[Wed Feb 15 17:19:41 2017].462013 msConnPoolClose(__rgealti_5m_courbes_niveau_CONTOUR__,0x7f45fda2ca00)
[Wed Feb 15 17:19:41 2017].462023 Entering msContourLayerReadRaster().
[Wed Feb 15 17:19:41 2017].462029 msContourLayerReadRaster(): Entering transform.
[Wed Feb 15 17:19:41 2017].462034 msContourLayerReadRaster(): src=138,0,665,123, dst=0,0,665,123
[Wed Feb 15 17:19:41 2017].467237 msConnPoolRegister(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7f45fda362e0)
[Wed Feb 15 17:19:41 2017].467257 msConnPoolRequest(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__) -> got 0x7f45fda362e0
[Wed Feb 15 17:19:41 2017].467404 msOGRFileWhichShapes: Setting spatial filter to 875716.95911215 6509443 878985.04088785 6511710
[Wed Feb 15 17:19:41 2017].467421 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].467440 msOGRFileNextShape: Returning shape=0, tile=0
[Wed Feb 15 17:19:41 2017].508823 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].508852 msOGRFileNextShape: Returning shape=1, tile=0
[Wed Feb 15 17:19:41 2017].508928 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].508944 msOGRFileNextShape: Returning shape=2, tile=0
[Wed Feb 15 17:19:41 2017].509051 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].509067 msOGRFileNextShape: Returning shape=3, tile=0
[Wed Feb 15 17:19:41 2017].509178 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].509199 msOGRFileNextShape: Returning shape=4, tile=0
[Wed Feb 15 17:19:41 2017].509661 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].509680 msOGRFileNextShape: Returning shape=5, tile=0
[Wed Feb 15 17:19:41 2017].509934 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].509953 msOGRFileNextShape: Returning shape=6, tile=0
[Wed Feb 15 17:19:41 2017].510161 Entering msContourLayerNextShape().
[Wed Feb 15 17:19:41 2017].510170 msOGRFileNextShape: Returning MS_DONE (no more shapes)
[Wed Feb 15 17:19:41 2017].510177 Entering msContourLayerClose().

So it considers 6 shapes - and here's the same log with the same shp2img query, but this time using DATA and pointing at the vrt:

[Wed Feb 15 17:22:56 2017].101479 Entering msContourLayerReadRaster().
[Wed Feb 15 17:22:56 2017].101519 msContourLayerReadRaster(): Entering transform.
[Wed Feb 15 17:22:56 2017].101526 msContourLayerReadRaster(): src=49138,24653,665,465, dst=0,0,665,465
[Wed Feb 15 17:22:56 2017].257137 msConnPoolRegister(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7fbc452edcb0)
[Wed Feb 15 17:22:56 2017].257360 msConnPoolRequest(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__) -> got 0x7fbc452edcb0
[Wed Feb 15 17:22:56 2017].257410 Entering msContourLayerWhichShapes().
[Wed Feb 15 17:22:56 2017].257416 msConnPoolRelease(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7fbc452edcb0)
[Wed Feb 15 17:22:56 2017].257421 msOGRLayerClose(__rgealti_5m_courbes_niveau_CONTOUR__).
[Wed Feb 15 17:22:56 2017].257425 msOGRFileClose(,0).
[Wed Feb 15 17:22:56 2017].257429 msConnPoolRelease(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7fbc452edcb0)
[Wed Feb 15 17:22:56 2017].257433 msConnPoolClose(__rgealti_5m_courbes_niveau_CONTOUR__,0x7fbc452edcb0)
[Wed Feb 15 17:22:56 2017].257446 Entering msContourLayerReadRaster().
[Wed Feb 15 17:22:56 2017].257454 msContourLayerReadRaster(): Entering transform.
[Wed Feb 15 17:22:56 2017].257459 msContourLayerReadRaster(): src=49138,24653,665,465, dst=0,0,665,465
[Wed Feb 15 17:22:56 2017].284745 msConnPoolRegister(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__,0x7fbc45305ff0)
[Wed Feb 15 17:22:56 2017].284770 msConnPoolRequest(rgealti_5m_courbes_niveau,__rgealti_5m_courbes_niveau_CONTOUR__) -> got 0x7fbc45305ff0
[Wed Feb 15 17:22:56 2017].284978 msOGRFileWhichShapes: Setting spatial filter to 875716.95911215 6509443 878985.04088785 6511710
[Wed Feb 15 17:22:56 2017].284996 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].285021 msOGRFileNextShape: Returning shape=0, tile=0
[Wed Feb 15 17:22:56 2017].285495 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].285514 msOGRFileNextShape: Returning shape=1, tile=0
[Wed Feb 15 17:22:56 2017].285556 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].285579 msOGRFileNextShape: Returning shape=2, tile=0
[Wed Feb 15 17:22:56 2017].286385 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].286403 msOGRFileNextShape: Returning shape=3, tile=0
[Wed Feb 15 17:22:56 2017].286474 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].286488 msOGRFileNextShape: Returning shape=4, tile=0
[Wed Feb 15 17:22:56 2017].286554 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].286570 msOGRFileNextShape: Returning shape=5, tile=0
[Wed Feb 15 17:22:56 2017].286755 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].286773 msOGRFileNextShape: Returning shape=6, tile=0
[Wed Feb 15 17:22:56 2017].287003 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].287031 msOGRFileNextShape: Returning shape=7, tile=0
[Wed Feb 15 17:22:56 2017].288078 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].288113 msOGRFileNextShape: Returning shape=8, tile=0
[Wed Feb 15 17:22:56 2017].289392 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].289420 msOGRFileNextShape: Returning shape=9, tile=0
[Wed Feb 15 17:22:56 2017].290358 Entering msContourLayerNextShape().
[Wed Feb 15 17:22:56 2017].290367 msOGRFileNextShape: Returning MS_DONE (no more shapes)
[Wed Feb 15 17:22:56 2017].290375 Entering msContourLayerClose().

Something seems wrong in the src/dst debug msg coming from https://github.com/mapserver/mapserver/blob/branch-7-0/mapcontour.c#L383 .. maybe a hint ?

landryb avatar Feb 15 '17 16:02 landryb

I'm experiencing the same problem in 7.0.6. Is there a chance this could get fixed?

ChadDean avatar Sep 08 '17 13:09 ChadDean

Anyone one have a test case? I have no access to tiled contour data... --Steve

sdlime avatar Sep 08 '17 20:09 sdlime

What kind of testcase do you need ? a wms server with a 'broken' contour layer ? Otherwise it's easy to setup a local server with open DEM data...

landryb avatar Sep 11 '17 07:09 landryb

Attached is a test case with a tileindex shapefile using the 2 DEMs below.

ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/NED/13/IMG/USGS_NED_13_n48w100_IMG.zip ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/NED/13/IMG/n49w100.zip

There are 2 extents in the mapfile. 1 extent attempting to show contours over both images and another outside the tile index that crashes Mapserver. ContourErrors.zip

ChadDean avatar Sep 11 '17 14:09 ChadDean

Looking at using contours again and it appears this still an issue. Any chance someone could take a look at it?

ChadDean avatar May 23 '18 13:05 ChadDean

I encounter the same issue on 7.6.4

pelord avatar Mar 23 '22 15:03 pelord

Just ran into this myself in version 8.0.1.

geographika avatar Sep 29 '23 07:09 geographika

Notes from investigations:

msDrawRasterIterateTileIndex only ever returns one file image at a time here.

In order for contours to be generated for multiple images in a tileindex the code would need to be refactored to use a loop similar to when creating raster layers, here

  done = MS_FALSE;
  while(done != MS_TRUE) {

    if(layer->tileindex) {
      status = msDrawRasterIterateTileIndex(layer, tlp, &tshp,
                                            tileitemindex, tilesrsindex,
                                            tilename, sizeof(tilename),
                                            tilesrsname, sizeof(tilesrsname));

Support for tileindexes was added after the initial contour feature was added via https://github.com/MapServer/MapServer/commit/946c210c2976dde1a2d10152655feaacedfda5cf "Contour layer: support tileindex and WMS time".

I'm unsure if this is by design - maybe support for tileindexes should only get the first tile in an extent, although at any image boundary there will be an issue.

geographika avatar Sep 29 '23 08:09 geographika

@rouault - unlikely you have any recollection of a commit from 7 years ago, but was the tileindex support in https://github.com/MapServer/MapServer/commit/946c210c2976dde1a2d10152655feaacedfda5cf just to stop it crashing in a tileindex layer and return a result for a single raster? I can't see an easy fix to get it working for multiple rasters in a tileindex - it looks like a major refactor of mapcountour.c would be needed.

geographika avatar Oct 04 '23 21:10 geographika

but was the tileindex support in 946c210 just to stop it crashing in a tileindex layer and return a result for a single raster?

no, but it was limited to the case of WMS TIME= support where the tileindex contains raster of the same extent but at different time positions. To contour a mosaic of rasters, the files of the tileindex should be put as sources of a single GDAL VRT, which would be the GDAL raster subject to contouring.

rouault avatar Oct 04 '23 22:10 rouault

Thanks that makes sense. I'll work up an example of a GDAL VRT.

geographika avatar Oct 04 '23 22:10 geographika

For reference the following approach works perfectly:

If you are creating a TILEINDEX using the following command:

gdaltindex index.shp myfolder/*.tif

With a LAYER as follows:

LAYER
    NAME "Contours"
    TYPE LINE
    CONNECTIONTYPE CONTOUR
    TILEINDEX "index.shp" # only one tile is loaded at a time
    TILEITEM "location"

    PROCESSING "BANDS=1"
    PROCESSING "CONTOUR_ITEM=elevation"
    PROCESSING "CONTOUR_INTERVAL=20"

You can replace it by creating a VRT file:

gdalbuildvrt output.vrt myfolder/*.tif

And the following LAYER definition:

LAYER
    NAME "Contours"
    TYPE LINE
    CONNECTIONTYPE CONTOUR
    DATA "output.vrt"

    PROCESSING "BANDS=1"
    PROCESSING "CONTOUR_ITEM=elevation"
    PROCESSING "CONTOUR_INTERVAL=20"

I think this workaround and the complications in trying to rewrite the contour layer code to support tileindexes will probably mean this issue is as fixed as it can be.

geographika avatar Oct 09 '23 12:10 geographika

i would agree with you since that's what i ended up doing to workaround this problem (see my first comment).. that's just a bit awkward when your tileindexes are in postgis (eg separate data from sidecar/tileindex files..).

landryb avatar Oct 09 '23 12:10 landryb