search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.87k stars 2.54k forks source link

Bad output resolution for GoogleMapsCompatible COG with space projections #5532

Closed Plantain closed 2 years ago

Plantain commented 2 years ago

Expected behavior and actual behavior.

Observe the significant loss in fidelity between left and right.

ex

Zoomed further: ex2

Left is the XYZ tile 7/116/76 as generated with gdalwarp -t_srs EPSG:3785 -tr 700 700 -r lanczos /gdal_translate to MBTILES/gdaladdo

Right is the same tile from the GoogleMapsCompatible COG generated via gdal_translate -of COG himawari-sat.tiff -co TILING_SCHEME=GoogleMapsCompatible -r lanczos -co BLOCKSIZE=512 -co ZOOM_LEVEL_STRATEGY=UPPER

The goal is to reproject GeoTIFF's of satellite imagery into GoogleMapsCompatible COG's, preserving the input resolution/fidelity

I suspect that GDALSuggestedWarpOutput is suggesting a bad output resolution because it is struggling with the space projection, which is not overridable when generating COG's like it is with most other outputs. The calculated pixelsize of 1222.99 is considerably less than the input resolution.

It should be possible to set a target zoom level, or specify an output resolution when generating a GoogleMapsCompatible COG, so that the correct resolution can be used on satellite data. Or, the correct output resolution should be calculated.

Original data:

gdalinfo himawari-sat.tiff 
Driver: GTiff/GeoTIFF
Files: himawari-sat.tiff
Size is 18408, 19008
Coordinate System is:
PROJCRS["unknown",
    BASEGEOGCRS["GCS_unknown",
        DATUM["D_unknown",
            ELLIPSOID["unknown",6378137,298.257024882273,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Geostationary Satellite (Sweep Y)"],
        PARAMETER["Longitude of natural origin",140.7,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Satellite Height",35785863,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-4701999.972949326038361,4501999.974099928513169)
Pixel Size = (499.999997123492847,-499.999997123492619)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2022:03:27 21:00:00
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-4701999.973, 4501999.974) 
Lower Left  (-4701999.973,-5001999.971) 
Upper Right ( 4501999.974, 4501999.974) 
Lower Right ( 4501999.974,-5001999.971) 
Center      (  -99999.999, -249999.999) (139d48' 2.90"E,  2d15'42.86"S)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha

Left (gdalwarp -t_srs EPSG:3785 -tr 700 700 -r lanczos):

gdalinfo himawari-3857-700.tiff
Driver: GTiff/GeoTIFF
Files: himawari-3857-700.tiff
Size is 10973, 10807
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        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 = (12245143.987260092049837,1118889.974857959663495)
Pixel Size = (700.000000000000000,-700.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2022:03:27 21:00:00
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (12245143.987, 1118889.975) (110d 0' 0.00"E, 10d 0' 0.00"N)
Lower Left  (12245143.987,-6446010.025) (110d 0' 0.00"E, 49d59'54.47"S)
Upper Right (19926243.987, 1118889.975) (179d 0' 1.78"E, 10d 0' 0.00"N)
Lower Right (19926243.987,-6446010.025) (179d 0' 1.78"E, 49d59'54.47"S)
Center      (16085693.987,-2663560.025) (144d30' 0.89"E, 23d15'38.03"S)
Band 1 Block=10973x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=10973x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=10973x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=10973x1 Type=Byte, ColorInterp=Alpha

Right:

gdalinfo himawari-cog.tiff
Driver: GTiff/GeoTIFF
Files: himawari-cog.tiff
Size is 32768, 14336
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        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 = (-20037508.342789243906736,7514065.628545965999365)
Pixel Size = (1222.992452562820063,-1222.992452562819835)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2022:03:27 21:00:00
Image Structure Metadata:
  INTERLEAVE=PIXEL
  LAYOUT=COG
Tiling Scheme:
  NAME=GoogleMapsCompatible
  ZOOM_LEVEL=6
Corner Coordinates:
Upper Left  (-20037508.343, 7514065.629) (180d 0' 0.00"W, 55d46'35.66"N)
Lower Left  (-20037508.343,-10018754.171) (180d 0' 0.00"W, 66d30'47.74"S)
Upper Right (20037508.343, 7514065.629) (180d 0' 0.00"E, 55d46'35.66"N)
Lower Right (20037508.343,-10018754.171) (180d 0' 0.00"E, 66d30'47.74"S)
Center      (       0.000,-1252344.271) (  0d 0' 0.01"E, 11d10'42.25"S)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224
  Mask Flags: PER_DATASET ALPHA 
  Overviews of mask band: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224
Band 4 Block=512x512 Type=Byte, ColorInterp=Alpha
  Overviews: 16384x7168, 8192x3584, 4096x1792, 2048x896, 1024x448, 512x224

Operating system

Ubuntu 21.10

GDAL version and provenance

This example was generated GDAL 3.2.2, released 2021/03/05 Same behaviour also observed with GDAL 3.4.1.

jratike80 commented 2 years ago

What resolution gdalwarp computes if you do not define it with -tr?

I tried it myself with a test image that I created with gdal_create gdal_create -a_srs space.prj -outsize 18408 19008 -a_ullr -4701999.973 4501999.974 4501999.974 -5001999.971 space.tif

gdalwarp -t_srs epsg:3857 -r lanczos space.tif space_warped.tif Gdalwarp prints errors ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds. but it does create an output that has a pixel size of Pixel Size = (1650.183828510722833,-1650.183828510722833)

Plantain commented 2 years ago

@jratike80 1650.183862530689794

I uploaded the original input data here if anyone wants to try: http://static2.skysight.io/demo/himawari-sat.tiff

gdalinfo out.tiff
Driver: GTiff/GeoTIFF
Files: out.tiff
Size is 24285, 10507
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        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 = (-20037508.325681626796722,7446961.240155051462352)
Pixel Size = (1650.183862530689794,-1650.183862530689794)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2022:03:27 21:00:00
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-20037508.326, 7446961.240) (180d 0' 0.00"W, 55d26' 9.83"N)
Lower Left  (-20037508.326,-9891520.603) (180d 0' 0.00"W, 66d 3'12.82"S)
Upper Right (20037206.776, 7446961.240) (179d59'50.25"E, 55d26' 9.83"N)
Lower Right (20037206.776,-9891520.603) (179d59'50.25"E, 66d 3'12.82"S)
Center      (    -150.775,-1222279.682) (  0d 0' 4.88"W, 10d54'47.99"S)
Band 1 Block=24285x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=24285x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=24285x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=24285x1 Type=Byte, ColorInterp=Alpha
Plantain commented 2 years ago

I also get similar errors: ERROR 1: Point outside of projection domain ERROR 1: Point outside of projection domain ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds. Warning 1: Unable to compute source region for output window 7156,0,7156,12384, skipping. ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds. Warning 1: Unable to compute source region for output window 7156,12384,7156,12385, skipping. ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds. Warning 1: Unable to compute source region for output window 14312,0,14313,24769, skipping. ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds. Warning 1: Unable to compute source region for output window 28625,0,7156,12384, skipping. ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds.

But for a COG, there is no way to manually specify output bounds/resolution with GoogleMapsCompatible.

vincentsarago commented 2 years ago

@Plantain maybe what you can do, is first create a VRT (with the resolution you want) and then create a COG 🤷(I know this doesn't fix the issue, if one, but it should give you want you want)

Plantain commented 2 years ago

@vincentsarago that's approximately what I do at the moment to work around this - gdal_warp to epsg:3857, then gdal_translate to a GoogleMapsCompatible COG. I don't completely understand the handling of the VRT as to whether this is actually introducing another reprojection step with loss of precision / wasted computation or not compared to how this 'should'/could work all in one step?

Incidentally I've found if I actually do the projection to epsg:3857 (gdalwarp -t_srs EPSG:3785 -tr 700 700 in.tiff out.tiff; gdal_translate -of COG out.tiff out_cog.tiff) rather than outputting a VRT and translating that, that the whole operation completes much faster, and I don't really understand why that would be because as I understand it that will definitely be doing two reprojections.

jratike80 commented 2 years ago

Web Mercator tiling scheme is using EPSG:3857 so why would there be two reprojections?

I believe that the real difficulty with this image is somehow related to the coordinate system. The size of the original image is 18408, 19008 and when it is warped to EPSG:3857 with defaults the size of the resulting image is 24285, 10507. If the width goes down from 19008 to 10507 then there must happen quality loss. I think this in not normal for gdalwarp.

The ZOOM_LEVEL_STRATEGY=UPPER already forces GDAL to do oversampling compared to what it considers to be the optimal resolution. Making an option to force even more oversampling feels like trying to fix an error that happens in some other place but maybe something like ZOOM_LEVEL_STRATEGY=UPPER+1 could be supported even in normal cases it only balloons the file size without any improvement to quality.

VRT idea is good. These commands give a COG with pixel size of 611.5 m (means the next Web Mercator zoom level).

gdalwarp -of VRT -tr 700 700 -t_srs epsg:3857 space.tif space.vrt                  
gdal_translate -of COG -co TILING_SCHEME=GoogleMapsCompatible -r lanczos -co BLOCKSIZE=512 -co ZOOM_LEVEL_STRATEGY=UPPER space.vrt space_cog4.tif
Plantain commented 2 years ago

Web Mercator tiling scheme is using EPSG:3857 so why would there be two projections?

Is GDAL clever enough pad the alignment for GoogleMapsCompatible (to make tiles align with Web Mercator tiles) without any more warping? I don't know enough about the internals of GDAL.

The ZOOM_LEVEL_STRATEGY=UPPER already forces GDAL to do oversampling compared to what it considers to be the optimal resolution. Making an option to force even more oversampling feels like trying to fix an error that happens in some other place

I think the root problem is GDAL doesn't know how to calculate the optimal resolution of space projected imagery? I agree that twiddling ZOOM_LEVEL_STRATEGY sounds like a hacky fix.

rouault commented 2 years ago

I think the root problem is GDAL doesn't know how to calculate the optimal resolution of space projected imagery?

that's a bit deeper than that. Neither the Geostationary Satellite projection nor the WebMercator projection are equal-area projection. Consequently the actual area on the ellipsoid of let's say a square of 700x700 meter is generally not 700x700 squared meters, except at the lines of true scale of the projections, and elsewhere it might be much larger or smaller. Consequently, and this is the general case even for equal-area projections, GDAL uses a strategy that tries to preserve the number of pixels, while outputing square pixels in the target projection : 18408 19008 = 349 899 264 and 32768 14336 = 469 762 048 are not so far in this regard. This is a compromise between getting the best quality possible (that could lead to an enormous output image) and preserving the average amount of information. The Geostationary Satellite projection, when used with images that incorporate portions of the sky, is of course a bad case for this heuristics as a significant portion of the pixels of the source image will not participate to the target image... That said being able to specify a desired target zoom level is a reasonable way forward.