3

Adding in Leaflet a WMTS layer with the 'L.tileLayer' to a COG Tiff results in an offset. Both services are in the same EPSG:2056 system

This code

<html>
<head>
    !-- Version 2.0 --> <!-- full stretch -->
    <meta name="viewport" content="width=device-width, initial-scale=1">
    <link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css"/>
    <style>
        #map {
            bottom: 0;
            left: 0;
            position: absolute;
            right: 0;
            top: 0;
        }
    </style>
</head>
<body>
<div id="map"></div>
<script src="https://unpkg.com/[email protected]/dist/leaflet.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.6.0/proj4.js"></script>
<script src="https://unpkg.com/georaster"></script>
<script src="https://unpkg.com/georaster-layer-for-leaflet"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4leaflet/1.0.2/proj4leaflet.min.js"></script>


<script>
    // Define custom CRS for EPSG 2056
    var crs2056 = new L.Proj.CRS('EPSG:2056',
        "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs",
        {
            resolutions: [1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 0.5, 0.25, 0.125, 0.0625],
            origin: [2485071.58, 1299941.79],
        }
    );

    // Initialize Leaflet map with custom CRS
    var map = L.map('map', {
        crs: crs2056,
        maxZoom: crs2056.options.resolutions.length,

    }).setView([46, 8], 1);

    // Add TMS Layer
    L.tileLayer('https://wmts.asit-asso.ch/wmts/1.0.0/asitvd.fond_pourortho/default/default/0/2056/{z}/{y}/{x}.png', {
    maxZoom: crs2056.options.resolutions.length,
    tileSize: 128
    }).addTo(map);
 
    // Add WMTS Layer: Does not work since  Tilematrix is wrongly set with 20
    //var map4 = new L.TileLayer.WMTS("https://wmts.asit-asso.ch/wmts", {
    //  layer: 'asitvd.fond_pourortho',
    //  style: "default",
    //  format: 'image/png',
    //  tilematrixSet: "2056"
        
    //}).addTo(map);

    const url_to_geotiff_file = new URLSearchParams(location.search).get("url");
    console.log("url_to_geotiff_file:", url_to_geotiff_file);

    if (!url_to_geotiff_file) {
        setTimeout(function() {
            // If no URL is provided, redirect to a sample GeoTIFF
            const parser = new URL(window.location);
            parser.searchParams.set("url", "https://d29cp2gnktw6by.cloudfront.net/data/S2_LEVEL_2A/20240101T102431/S2_L2A_SR_20240101T102431_20240101T102428_10M_run20240102.tif");
            window.location = parser.href;
        }, 2 * 1000);
    } else {
        console.log("url_to_geotiff_file:", url_to_geotiff_file);
        // Parse the GeoTIFF data
        parseGeoraster(url_to_geotiff_file).then(function (georaster) {
            console.log("georaster:", georaster);
            // based on https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/l2a_optimized/    
        
        // Contrast enhance / highlight compress
        const maxR = 3.0; // max reflectance
        const midR = 0.13;
        const sat = 1.2;
        const gamma = 1.8;
        
        
        function evaluatePixel(val0,val1,val2) {
            //console.log("val0:", val0, "val1:", val1, "val2:", val2);
            const rgbLin = satEnh(sAdj(val0), sAdj(val1), sAdj(val2));
            //console.log("rgbLin:", rgbLin);
            //console.log("sRGB0:", sRGB(rgbLin[0]), "sRGB1:", sRGB(rgbLin[1]), "sRGB2:", sRGB(rgbLin[2]));
            return [sRGB(rgbLin[0]), sRGB(rgbLin[1]), sRGB(rgbLin[2])];
        }
        
        function sAdj(a) {
            return adjGamma(adj(a, midR, 1, maxR));
        }

        
        const gOff = 0.01;
        const gOffPow = Math.pow(gOff, gamma);
        const gOffRange = Math.pow(1 + gOff, gamma) - gOffPow;


        function adjGamma(b) {
            return (Math.pow((b + gOff), gamma) - gOffPow) / gOffRange;
        }
        
        // Saturation enhancement
        function satEnh(r, g, b) {
            const avgS = (r + g + b) / 3.0 * (1 - sat);
            return [clip(avgS + r * sat), clip(avgS + g * sat), clip(avgS + b * sat)];
        }

        function clip(s) {
            return s < 0 ? 0 : s > 1 ? 1 : s;
        }

        //contrast enhancement with highlight compression
        function adj(a, tx, ty, maxC) {
            var ar = clip(a / maxC, 0, 1);
            return ar * (ar * (tx / maxC + ty - 1) - ty) / (ar * (2 * tx / maxC - 1) - tx / maxC);
        }

        const sRGB = (c) => c <= 0.0031308 ? (12.92 * c) : (1.055 * Math.pow(c, 0.41666666666) - 0.055);


            const layer = new GeoRasterLayer({
                attribution: "Contains modified Copernicus Sentinel data 2024",
                debugLevel: 0,
                georaster,
                resolution: 128,
                opacity:1,
                resampleMethod: "nearest",
                pixelValuesToColorFn: values => {
                    // Check if any of the bands have a value of 9999 or 0
                    if (values.some(value => value === 9999 )) {
                    //if (values.some(value => value === 9999 || value === 0)) {
                        // If any band has a value of 9999 or 0, return transparent color
                        return 'rgba(0, 0, 0, 0)'; // Transparent color
                    } else {
                        // Apply min-max normalization to each band
                        const r = Math.round(evaluatePixel(values[0]/10000, values[1]/10000, values[2]/10000)[0] * 255);
                        const g = Math.round(evaluatePixel(values[0]/10000, values[1]/10000, values[2]/10000)[1] * 255);
                        const b = Math.round(evaluatePixel(values[0]/10000, values[1]/10000, values[2]/10000)[2] * 255);


                        // Return RGB color
                        return `rgba(${r},${g},${b},1)`;
                    }
                }
            });

            layer.addTo(map);
            map.fitBounds(layer.getBounds());
        });
    }
</script>
</body>
</html>

results in this map

enter image description here

any ideas why the 'asitvd.fond_pourortho' tile layer is off in the up left corner?

1 Answer 1

5

The main problem is in the tile grid. WMTS layer from https://wmts.asit-asso.ch/ has the same projection, but different tile grid. There are some maps available from asit-asso.ch on https://asit-asso.ch/service/fonds-de-plan. Map is in OpenLayers. Source code inspection in browse debugger revealed the following about tile grid:

const ASITResolutions = [4000.0, 3750.0, 3500.0, 3250.0, 3000.0, 2750.0, 2500.0, 2250.0, 2000.0, 1750.0, 1500.0, 1250.0, 1000.0, 750.0, 650.0, 500.0, 250.0, 100.0, 50.0, 20.0, 10.0, 5.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.25, 0.1, 0.05];
const ASITMatrixIds = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29];
const WMTSOrigin = [2420000.0, 1350000.0];

These values have to be used on projection definition.

One other thing is that WMTS tile size is standard 256, so tileSize: 128 option should be left out.

Relevant part of the code then looks something like this:

var crs2056 = new L.Proj.CRS('EPSG:2056',
   "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs",
   {
     resolutions: [4000.0, 3750.0, 3500.0, 3250.0, 3000.0, 2750.0, 2500.0, 2250.0, 2000.0, 1750.0, 1500.0, 1250.0, 1000.0, 750.0, 650.0, 500.0, 250.0, 100.0, 50.0, 20.0, 10.0, 5.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.25, 0.1, 0.05],
     origin: [2420000.0, 1350000.0]
 });

var map = L.map('map', {
  crs: crs2056
}).setView([46, 8], 1);

L.tileLayer('https://wmts.asit-asso.ch/wmts/1.0.0/asitvd.fond_pourortho/default/default/0/2056/{z}/{y}/{x}.png', {
  maxZoom: crs2056.options.resolutions.length,
  zIndex: 5
}).addTo(map);

This is then the result:

enter image description here

1
  • I was looking at this service in QGIS, I can can get the KVP service to work, but all the ZYX requests are invalid, though hand crafted comparable requests return results. I guess that QGIS is assuming one grid set for the zoom level and sending inappropriate requests.
    – nmtoken
    Commented Feb 23 at 13:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.