Opened 6 years ago
Closed 6 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
$ 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 , 6 years ago
comment:2 by , 6 years ago
One way to improve things too is to add "-r bilinear" in the gdalbuildvrt invokation
comment:3 by , 6 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?
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 , 6 years ago
| Resolution: | → invalid |
|---|---|
| Status: | new → closed |
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
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:
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.