GeoTIFF reproject

[panterlo]

I have testing out the features of GeoTIFF and have a GeoTIFF COG with the following data:

PROJCRS["SWEREF99 TM",
    BASEGEOGCRS["SWEREF99",
        DATUM["SWEREF99",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4619]],
    CONVERSION["SWEREF99 TM",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",15,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (medium and small scale)."],
        AREA["Sweden - onshore and offshore."],
        BBOX[54.96,10.03,69.07,24.17]],
    ID["EPSG",3006]]
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 7813x7813, 3907x3907, 1954x1954, 977x977, 489x489
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 7813x7813, 3907x3907, 1954x1954, 977x977, 489x489
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 7813x7813, 3907x3907, 1954x1954, 977x977, 489x489
Band 4 Block=512x512 Type=Byte, ColorInterp=Undefined
  Overviews: 7813x7813, 3907x3907, 1954x1954, 977x977, 489x489

Putting this into Martin gives this result:

Running:

gdalwarp -t_srs EPSG:3857 -of cog -dstalpha -srcnodata 0 -ot BYTE -b 1 -b 2 -b 3 -overwrite -multi --config GDAL_CACHEMAX 9000 -wm 9000 -co BLOCKSIZE=256 -co BIGTIFF=YES -co COMPRESS=DEFLATE -co PREDICTOR=2 -co LEVEL=9 -co NUM_THREADS=ALL_CPUS test.tif reprojected.tif

Gives the following result where the tile is repeated and not matched to the map:

Loading both GeoTIFF into QGIS all looks correct. So can someone help out here to understand how the GeoTIFF should look like when they are loaded into Martin ?

[sharkAndshark]

Hi @panterlo I will look into this later would you mind share the file for debugging?

[panterlo]

Of course, here in a GeoTIFF you can download (this is a file with four bands and I haven't run any gdalwarp command on this file):

https://drive.google.com/drive/folders/12XRRaMceXZXCJRdrfJxrbdMu97kr_Ubv?usp=sharing

[sharkAndshark]

Hi @panterlo
I reprojected this and get it served by martin with a maptiler basemap.

gdalwarp -t_srs EPSG:3857 -of COG ./16_o72775_6950_25_mr24.tif  reprojected.tif

Driver: GTiff/GeoTIFF
Files: ./reprojected.tif
Size is 16634, 16654
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (2140275.216185077559203,9763663.895291147753596)
Pixel Size = (0.386031731321650,-0.386031731321650)
Metadata:
  AREA_OR_POINT=Area
  OVR_RESAMPLING_ALG=NEAREST
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
  LAYOUT=COG
Corner Coordinates:
Upper Left  ( 2140275.216, 9763663.895) ( 19d13'35.11"E, 65d34'59.13"N)
Lower Left  ( 2140275.216, 9757234.923) ( 19d13'35.11"E, 65d33'33.15"N)
Upper Right ( 2146696.468, 9763663.895) ( 19d17' 2.77"E, 65d34'59.13"N)
Lower Right ( 2146696.468, 9757234.923) ( 19d17' 2.77"E, 65d33'33.15"N)
Center      ( 2143485.842, 9760449.409) ( 19d15'18.94"E, 65d34'16.15"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 8318x8328, 4159x4164, 2080x2083, 1040x1041, 521x521
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 8318x8328, 4159x4164, 2080x2083, 1040x1041, 521x521
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 8318x8328, 4159x4164, 2080x2083, 1040x1041, 521x521
Band 4 Block=512x512 Type=Byte, ColorInterp=Undefined
  Overviews: 8318x8328, 4159x4164, 2080x2083, 1040x1041, 521x521

import './style.css';
import { Map, View } from 'ol';
import TileLayer from 'ol/layer/Tile';
import OSM from 'ol/source/OSM';
import XYZ from 'ol/source/XYZ.js';


import TileGrid from 'ol/tilegrid/TileGrid.js';

var custom_grid = new TileGrid({
  extent: [
    2140275.2161850776,9757234.922837717,2146696.468003882,9763663.895291148
    ],
  resolutions: [
    12.32485953705246,6.174280595004165,3.0871402975020823,1.5439412884838497,0.7719706442419249,0.3860317313216503
 ],
  tileSize: [512, 512],
  origin: [2140275.2161850776,9763663.895291148]
});

var source = new XYZ({
  url: "http://localhost:3000/reprojected/{z}/{x}/{y}",
  tileGrid: custom_grid
});

const map = new Map({
  target: 'map',
  layers: [
    new TileLayer({
      source: new XYZ({
        url: 'https://api.maptiler.com/maps/basic/{z}/{x}/{y}.png?key=vgdy5zgzM4IIqQCC6SOn',
        attributions: '© MapTiler © OpenStreetMap contributors'
      })
    }),
    new TileLayer({
      source: source
    }),

  ],
  view: new View({
    center: [2140275.2161850776,9763663.895291148],
    zoom: 14
  })
});

Seems the tiles are where they should be?

To load the tiles from COG, the custom tilegrid is required. But sadly we are still trying to figure out a better way to output them and not break the compatibility with TileJSON.

More discussion could be found in #343 either. I made a proposal there.

[sharkAndshark]

Before we know the best way to tell user the custom tilegrid of COG, at least maybe we could make a CLI like

cog world_cities.tif
File Path: /path/to/world_cities.tif
File Size: 10Mb
EPSG Code: 3857
Tile Grid: 
    origin: ....
    extent: ....
    resolutions: ....
    tile schema: ....

[sharkAndshark]

And another solution is allow tile clipping functionality in Martin to align COG (Cloud Optimized GeoTIFF) tiles with a user-specified grid system.

[panterlo]

Thanks !