Opened 5 years ago

Closed 5 years ago

#7164 closed defect (invalid)

vrt multiresolution mosaic glitch

Reported by: mclausen Owned by: warmerdam
Priority: normal Milestone:
Component: default Version: unspecified
Severity: normal Keywords:
Cc:

Description

we are trying to setup an aerial image multi resolution mosaic with vrt. These vrt files will be used afterwards for the visualization with mapserver 6.4. we would like to combine two raster datasets one with 10cm pixel size and one with 25cm pixel size in one vrt file. the aerial image data is in the same reference system.

what we observe is, if you overlay the mixed resolution vrt (dop25_dop10_kombo.vrt) with the lower resolution vrt (dop25.vrt) or with the source dataset you can see an offset in the lower resolution parts of the image. you can use this link to visualize the offset [1]. The offset is not evenly distributed on the image and it depends on zoomlevel and image detail. The same effect is visible in gqis 2.14.8 and arcmap 10.4.1 using this vrt files directly.

do you have any idea what is causing this effect and how we could get around it?

you will find the data here [2] this is how we have created the vrt files and some metadata about the source images [3]

thank you in advance for your help.

marcel

[1] https://s.geo.admin.ch/7719c8ee5e

[2] https://data.geo.admin.ch/test/ltclm/testdata.tar.gz

[3]

$ echo 1126-33.tif DOP10_LV95_25861219_2015_1_15.tif | xargs gdalbuildvrt -resolution user -tr 0.1 0.1 -srcnodata "0 0 0"  dop25_dop10_kombo.vrt 
$ gdalbuildvrt dop25.vrt 1126-33.tif 
$ gdalbuildvrt dop10.vrt DOP10_LV95_25861219_2015_1_15.tif  
$ gdalbuildvrt --version
GDAL 2.3.0dev, released 2017/99/99

source image information 10cm pixel size:

$ gdalinfo DOP10_LV95_25861219_2015_1_15.tif
Driver: GTiff/GeoTIFF
Files: DOP10_LV95_25861219_2015_1_15.tif
Size is 10000, 10000
Coordinate System is:
PROJCS["CH1903+ / LV95",
    GEOGCS["CH1903+",
        DATUM["CH1903+",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[674.374,15.056,405.346,0,0,0,0],
            AUTHORITY["EPSG","6150"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4150"]],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",46.95240555555556],
    PARAMETER["longitude_of_center",7.439583333333333],
    PARAMETER["azimuth",90],
    PARAMETER["rectified_grid_angle",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",2600000],
    PARAMETER["false_northing",1200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","2056"]]
Origin = (2586000.000000000000000,1220000.000000000000000)
Pixel Size = (0.100000000000000,-0.100000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2014:01:07 13:45:26
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS4 Windows
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 2586000.000, 1220000.000) (  7d15'18.11"E, 47d 7'55.84"N)
Lower Left  ( 2586000.000, 1219000.000) (  7d15'18.22"E, 47d 7'23.46"N)
Upper Right ( 2587000.000, 1220000.000) (  7d16' 5.57"E, 47d 7'55.92"N)
Lower Right ( 2587000.000, 1219000.000) (  7d16' 5.67"E, 47d 7'23.53"N)
Center      ( 2586500.000, 1219500.000) (  7d15'41.89"E, 47d 7'39.69"N)
Band 1 Block=10000x1 Type=Byte, ColorInterp=Red
Band 2 Block=10000x1 Type=Byte, ColorInterp=Green
Band 3 Block=10000x1 Type=Byte, ColorInterp=Blue

source image information 25cm pixel size

$ gdalinfo 1126-33.tif
Driver: GTiff/GeoTIFF
Files: 1126-33.tif
       1126-33.tif.ovr
Size is 17500, 12000
Coordinate System is:
PROJCS["CH1903+ / LV95",
    GEOGCS["CH1903+",
        DATUM["CH1903+",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[674.374,15.056,405.346,0,0,0,0],
            AUTHORITY["EPSG","6150"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4150"]],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",46.95240555555556],
    PARAMETER["longitude_of_center",7.439583333333333],
    PARAMETER["azimuth",90],
    PARAMETER["rectified_grid_angle",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",2600000],
    PARAMETER["false_northing",1200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","2056"]]
Origin = (2585000.000000000000000,1221000.000000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=3 (pixels/cm)
  TIFFTAG_XRESOLUTION=28.346458
  TIFFTAG_YRESOLUTION=28.346458
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  ( 2585000.000, 1221000.000) (  7d14'30.54"E, 47d 8'28.15"N)
Lower Left  ( 2585000.000, 1218000.000) (  7d14'30.90"E, 47d 6'51.00"N)
Upper Right ( 2589375.000, 1221000.000) (  7d17'58.19"E, 47d 8'28.46"N)
Lower Right ( 2589375.000, 1218000.000) (  7d17'58.45"E, 47d 6'51.30"N)
Center      ( 2587187.500, 1219500.000) (  7d16'14.52"E, 47d 7'39.74"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 8750x6000, 4375x3000, 2188x1500, 1094x750, 547x375, 274x188, 137x94
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 8750x6000, 4375x3000, 2188x1500, 1094x750, 547x375, 274x188, 137x94
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 8750x6000, 4375x3000, 2188x1500, 1094x750, 547x375, 274x188, 137x94

Change History (4)

comment:1 by Even Rouault, 5 years ago

I'm not sure it is an issue in GDAL itself. I believe the issue is more in the software using it, not using the GDALRasterIOEx() interface with subpixel accuracy and non-nearest resampling. I know that QGIS requests GDAL with integer source window coordinates, and with nearest resampling. Looking at QGIS requests to GDAL, there are also cases where it doesn't even request the same output size (e.g 812x838 on the original 25cm, and 811x838 on the 10cm VRT). All this combined is the likely cause for such effects.

I've experimented the following:

gdal_translate 1126-33.tif extract_ori.tif  -projwin 2587086.000 1219604.750 2587289.000 1219395.250
gdal_translate dop25_dop10_kombo.vrt extract_vrt_projwin.tif -projwin 2587086.000 1219604.750 2587289.000 1219395.250 -outsize 812 838 
gdal_translate dop25_dop10_kombo.vrt extract_vrt_projwin_bilinear.tif -projwin 2587086.000 1219604.750 2587289.000 1219395.250 -outsize 812 838 -r bilinear

So extracting the same area (exact on 25 cm resolution) with the same output size

If you display all them in QGIS at full resolution, you'll notice a shift between extract_ori.tif and extract_vrt_projwin.tif. But extract_ori.tif and extract_vrt_projwin_bilinear.tif will overlay fine. Which shows that if requested properly GDAL should give good results.

Another potential issue (which will show if you display at scales where 1 pixel > 25 cm) is that the overviews generated on 1126-33.tif probably use nearest resampling, which will cause subpixel shifts. If you generate with "gdaladdo -r bilinear", things will get better.

comment:2 by Even Rouault, 5 years ago

One way to improve things too is to add "-r bilinear" in the gdalbuildvrt invokation

comment:3 by mclausen, 5 years ago

thank you very much for the helpful hints. the overviews of the source raster data has been with gauss resampling

gdaladdo -ro -r gauss --config PHOTOMETRIC_OVERVIEW YCBCR --config COMPRESS_OVERVIEW JPEG --config JPEG_QUALITY_OVERVIEW 85 $i 2 4 8 16 32 64 128

using mapserver 6.4 we were able to get completely rid of this effect by using tileindex instead of vrt. Do you know what's the difference in the raster access between tileindex and vrt?

https://cdn.knightlab.com/libs/juxtapose/latest/embed/index.html?uid=d1db6a50-db1e-11e7-b263-0edaf8f81e27

the effect is also disappearing for vrt when using mapserver 7.0.4, 7.0.7 or 7.1dev. obviously mapserver has changed the gdal interface in the meantime.

best regards marcel

comment:4 by Even Rouault, 5 years ago

Resolution: invalid
Status: newclosed

I wouldn't trust gauss resampling to give sub-pixel accuracy.

Regarding tileindex vs vrt, with tileindex, MapServer manages itself the compositing of all images and can probably better handle the different resolutions rather than in the VRT which uses by default nearest resampling.

I'm not completely clear why MapServer 7.x would be better than MapServer 6.x regarding VRT access. The access to GDAL hasn't fundamentally changed, but oh well, if that works...

So closing that ticket

Note: See TracTickets for help on using tickets.