search

Deadwood-ai / deadwood-api

Main FastAPI application for the deadwood backend
GNU General Public License v3.0
0 stars 0 forks source link

rio_cogeo stuck writing cogs with 4 bands #40

Open JesJehle opened 2 months ago

JesJehle commented 2 months ago

In some cases, conversion of GeoTIFFs to COGs is stuck without any error during the writing of the COG. https://github.com/Deadwood-ai/deadwood-api/issues/29#issuecomment-2312839087

To reproduce:

converting it with gdal directly works:

gdal_translate debug-cog.tif debug-cog-result1.tif -of COG -co COMPRESS=JPEG -co QUALITY=75 -co OVERVIEWS=8

converting with reo_cogeo gets stuck: I'm using the docker environment here

reproduce the error

from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles

# Define the destination keyword arguments for COG creation
config = dict(
        # GDAL_NUM_THREADS="ALL_CPUS",
        GDAL_NUM_THREADS="2",
        GDAL_TIFF_INTERNAL_MASK=True,
        # GDAL_TIFF_OVR_BLOCKSIZE=f"{blocksize}",
    )

# Define COG profile
cog_profile = cog_profiles.get("jpeg")

# Call cog_translate with the required arguments
cog_translate(
    "/data/archive/bbac672f-749d-4644-a499-1d2021633819_debug-cog.tif",  # Source file
    "/data/cogs/test-process.tif",  # Destination file
    cog_profile,  # COG profile
    config, 
    overview_level=8,
    use_cog_driver=True # Destination keyword arguments
)

The output is:

/usr/local/lib/python3.12/site-packages/rio_cogeo/cogeo.py:241: UserWarning: Nodata/Alpha band will be translated to an internal mask band.
  warnings.warn(
Reading input: /data/archive/bbac672f-749d-4644-a499-1d2021633819_debug-cog.tif
  [####################################]  100%          
Adding overviews...
Updating dataset tags...
Writing output to: /data/cogs/test-process.tif

Writing does not make any progress.

ISSUE:

It seems to be because of the additional alpha band of the geotif, geotifs with only 3 bands work fine. I think this is due to JPEG compression, the rio-cogeo docs states that the JPEG compression profile can only be used for 3 band data.

FIX

The simplest fix would be to always ignore the 4th band and only use the first 3. But for some reason the gdal only approach works fine, so we would just have to switch back to the old implementation:

https://github.com/Deadwood-ai/deadwood-api/blob/8aec512dd535e33883541ddc45f501d1e16be13e/src/deadwood/cog.py#L22

@mmaelicke why did you use the rio implementation in the first place, any reason not to go back? What should we do?

cmosig commented 2 months ago

Thanks for debugging this :+1: Not sure about that library but usually everything is using gdal in the background anyway.

cmosig commented 2 months ago

Die lib nutzt rasterio im Hintergrund (also auch gdal). Vielleicht kann man die ersten drei Bänder der tif mit rasterio einlesen (rasterio.open(bands=[1,2,3])) und das object der lib übergeben.

Oder direkt in rasterio implementieren...spart man sich eine lib und einen point of failure

mmaelicke commented 2 months ago

Haben wir ja alles schon gemacht. Das Problem ist dass es es nicht ganz transparent ist welche flugs cogeo nutzt. Ohne die lib kann Janusch die cogs nicht ohne Fragmente und mit Alpha-Kanal darstellen. Ich wollte eigentlich die verschiedenen flags durchprobieren, aber die cogs mussten noch am vor dem Urlaub richtig dargestellt werden können, daher kam nur die schnellste Lösung in frage. Das Problem tauchte zum ersten Mal nach ca 150 uploads auf und ich musste gleichzeitig die prozessierung ändern ( queue etc), daher ging ich davon aus, dass die Probleme am neuen Code und nicht den Bildern liegen.

Ich bin immer dafür Probleme zu verstehen und nachhaltig zu lösen, vor allem eins nach dem anderen, das ist dann nur nicht immer auch der schnellste Weg.

Alles was ic testen konnte ist, dass nur jpeg Komprimierung Sinn macht. Alles andere ist zu groß oder kann nicht dargestellt werden (oder letztlich eine Variante, die auch jpeg verwendet)

mmaelicke commented 2 months ago

Dh wenn gdal funktioniert und du @JesJehle auch die Darstellung getestet hast, dann kannst du einfach den gdal Befehl im api Code updaten und wieder den gdal code verwenden. Was genau musste denn geändert werden?

JesJehle commented 2 months ago

Ich habe noch keine funktionierende Konfiguration der _rio_calculate_cog Funktion gefunden, mir ist nur aufgefallen, dass _gdal_calculate_cog direkt funktioniert. Allerdings habe ich mir die COGs noch nicht mit OpenLayers angeschaut.

Ich überprüfe jetzt also mal den Output der _gdal_calculate_cog mit GeoTIFs mit 3 und 4 bändern in OpenLayers bezüglich der no data values. Wenn das funktioniert, mache in ein PR um _gdal_calculate_cog zu verwenden.

mmaelicke commented 2 months ago

Genau das meinte ich. Die gdal funktioniert, die haben wir ja auch lange verwendet. Allerdings hat die Darstellung nicht geklappt. Mit systematischen Tests werden wir die richtige Config auch finden, wir müssen uns halt dafür die Zeit nehmen. Ich bin mir nicht sicher welches flag fehlt, bzw. falsch gesetzt ist

JesJehle commented 2 months ago

Also wenn dich das GeoTIFF mit 4 bändern mit der _gdal_calculate_cog Funktion prozessieren kann ich es mit OpenLayers visualisieren ohne Artefakte. https://codesandbox.io/p/sandbox/cog-forked-fkt4mp?from-embed=

image

Infos zu GeoTIFF:

Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
  Unit Type: metre
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
  Unit Type: metre
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
  Unit Type: metre
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Unit Type: metre

Infos zum processierten COG:

Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 23524x14195, 11762x7097, 5881x3548, 2940x1774, 1470x887, 735x443, 367x221, 183x110
  Mask Flags: PER_DATASET
  Overviews of mask band: 23524x14195, 11762x7097, 5881x3548, 2940x1774, 1470x887, 735x443, 367x221, 183x110
  Unit Type: metre
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 23524x14195, 11762x7097, 5881x3548, 2940x1774, 1470x887, 735x443, 367x221, 183x110
  Mask Flags: PER_DATASET
  Overviews of mask band: 23524x14195, 11762x7097, 5881x3548, 2940x1774, 1470x887, 735x443, 367x221, 183x110
  Unit Type: metre
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 23524x14195, 11762x7097, 5881x3548, 2940x1774, 1470x887, 735x443, 367x221, 183x110
  Mask Flags: PER_DATASET
  Overviews of mask band: 23524x14195, 11762x7097, 5881x3548, 2940x1774, 1470x887, 735x443, 367x221, 183x110
  Unit Type: metre

Die _gdal_calculate_cog scheint also einfach das 4. Band rauszunehmen. Ich kann allerdings auch keine information über no data values finden. Ich teste das mit unterschiedihcne Geotifs ma aus.

JesJehle commented 2 months ago

Ok, using a 3 band TIFF:

Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  NoData Value=255
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  NoData Value=255
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  NoData Value=255

conver to COG:

Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  NoData Value=255
  Overviews: 12666x12025, 6333x6012, 3166x3006, 1583x1503, 791x751, 395x375, 197x187, 98x93
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  NoData Value=255
  Overviews: 12666x12025, 6333x6012, 3166x3006, 1583x1503, 791x751, 395x375, 197x187, 98x93
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  NoData Value=255
  Overviews: 12666x12025, 6333x6012, 3166x3006, 1583x1503, 791x751, 395x375, 197x187, 98x93

It runs fine, but instead of a transparent background we all is white again. https://codesandbox.io/p/sandbox/cog-forked-rg2lwk?file=%2Fmain.js%3A19%2C1 image

Funny enough If i process the same TIFf with the `

i get the transparent background: image

The cog info is:

Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 12667x12025, 6334x6013, 3167x3007, 1584x1504, 792x752, 396x376, 198x188, 99x94
  Mask Flags: PER_DATASET
  Overviews of mask band: 12667x12025, 6334x6013, 3167x3007, 1584x1504, 792x752, 396x376, 198x188, 99x94
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 12667x12025, 6334x6013, 3167x3007, 1584x1504, 792x752, 396x376, 198x188, 99x94
  Mask Flags: PER_DATASET
  Overviews of mask band: 12667x12025, 6334x6013, 3167x3007, 1584x1504, 792x752, 396x376, 198x188, 99x94
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 12667x12025, 6334x6013, 3167x3007, 1584x1504, 792x752, 396x376, 198x188, 99x94
  Mask Flags: PER_DATASET
  Overviews of mask band: 12667x12025, 6334x6013, 3167x3007, 1584x1504, 792x752, 396x376, 198x188, 99x94

The difference seems to be the Mask Flags: PER_DATASET

mmaelicke commented 2 months ago

Exactly. Deswegen sind wir auf cogeo gegangen als quick n'dirty. Das war der Ausgangspunkt der ganzen Übung. Ich hab dann cogeo genommen obwohl es nur mit ein paar Bildern getestet wurde. Danach waren die Labels wichtiger.

Ich finde cogeo zwar hilfreich, aber schlecht dokumentiert. Es fehlt vermutlich nur eine Einstellung. Oder evtl kann man immer nur die ersten drei Bänder einlesen und als file-Like-object übergeben. Alternativ müssen wir mit gdal die cogeo Einstellungen nachbauen. Soweit ich weiß, setzt der den alphakanal als Maske für na werte ein, bin mir da aber nicht sicher.

Wenn ein Wert als na in den Metadaten steht muss man mit der jpeg Komprimierung aufpassen, weil das die pixelwerte ändert. Evtl. ist genau das das Problem. Evtl. sind die weißen Flächen eben nicht genau 0 und 0 aber der na wert (ist leider üblich bei cogs)

JesJehle commented 2 months ago

no completly shure but i found a configuration where _rio_calculate_cog works: the overview seems to be the problem.

    cog_translate(
        tiff_file_path,
        cog_target_path,
        output_profile,
        config=config,
        web_optimized=True,
        # overview_level=overviews, 
        indexes=(1, 2, 3),
        add_mask=True,
        use_cog_driver=True,
    )

I will test it with more TIFFs after the break.

COG looks fine image The 3 band image also seems to work: image no data values are transparent and the process runs smoothly

JesJehle commented 2 months ago

Ok, its only the manual set overview which cause the problem. This setting seems to work smoothly:

    cog_translate(
        tiff_file_path,
        cog_target_path,
        output_profile,
        config=config,
        web_optimized=True,
        # overview_level=overviews, 
        # indexes=(1, 2, 3),
        # add_mask=True,
        use_cog_driver=True,
    )

works for 4band images image and 3band images image

Docs about web_optimized:

    `web_optimized` now only serve to activate a default `WebMercatorQuad` tile matrix set.

Also, ich if dont spedify web_optimized the:

Driver: GTiff/GeoTIFF
Files: 869625ba-98a9-40e3-9a33-dc2de888df2b_debug-cog-3band_cog_jpeg_ovr8_q75.tif
Size is 25333, 24050
Coordinate System is:
GEOGCRS["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]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (8.623504895271152,48.024892172470260)
Pixel Size = (0.000000101680195,-0.000000101680195)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  JPEGTABLESMODE=1
  JPEG_QUALITY=75
  LAYOUT=COG
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (   8.6235049,  48.0248922) (  8d37'24.62"E, 48d 1'29.61"N)
Lower Left  (   8.6235049,  48.0224468) (  8d37'24.62"E, 48d 1'20.81"N)
Upper Right (   8.6260808,  48.0248922) (  8d37'33.89"E, 48d 1'29.61"N)
Lower Right (   8.6260808,  48.0224468) (  8d37'33.89"E, 48d 1'20.81"N)
Center      (   8.6247928,  48.0236695) (  8d37'29.25"E, 48d 1'25.21"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  NoData Value=255
  Overviews: 12666x12025, 6333x6012, 3166x3006, 1583x1503, 791x751, 395x375, 197x187, 98x93
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  NoData Value=255
  Overviews: 12666x12025, 6333x6012, 3166x3006, 1583x1503, 791x751, 395x375, 197x187, 98x93
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  NoData Value=255
  Overviews: 12666x12025, 6333x6012, 3166x3006, 1583x1503, 791x751, 395x375, 197x187, 98x93

If i specie it, however, i get the right projection, and it seems specifically for visualization.

Driver: GTiff/GeoTIFF
Files: 869625ba-98a9-40e3-9a33-dc2de888df2b_debug-cog-3band_cog_jpeg_ovr8_q75-web-optimized-only-3band.tif
Size is 15616, 22016
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 = (959963.083604980260134,6110997.091971814632416)
Pixel Size = (0.018661383858686,-0.018661383858686)
Metadata:
  AREA_OR_POINT=Area
  OVR_RESAMPLING_ALG=NEAREST
  TILING_SCHEME_NAME=WebMercatorQuad
  TILING_SCHEME_ZOOM_LEVEL=23
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  JPEGTABLESMODE=1
  JPEG_QUALITY=75
  LAYOUT=COG
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (  959963.084, 6110997.092) (  8d37'24.58"E, 48d 1'29.61"N)
Lower Left  (  959963.084, 6110586.243) (  8d37'24.58"E, 48d 1'20.73"N)
Upper Right (  960254.500, 6110997.092) (  8d37'34.01"E, 48d 1'29.61"N)
Lower Right (  960254.500, 6110586.243) (  8d37'34.01"E, 48d 1'20.73"N)
Center      (  960108.792, 6110791.667) (  8d37'29.29"E, 48d 1'25.17"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 7808x11008, 3904x5504, 1952x2752, 976x1376, 488x688
  Mask Flags: PER_DATASET
  Overviews of mask band: 7808x11008, 3904x5504, 1952x2752, 976x1376, 488x688
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 7808x11008, 3904x5504, 1952x2752, 976x1376, 488x688
  Mask Flags: PER_DATASET
  Overviews of mask band: 7808x11008, 3904x5504, 1952x2752, 976x1376, 488x688
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 7808x11008, 3904x5504, 1952x2752, 976x1376, 488x688
  Mask Flags: PER_DATASET
  Overviews of mask band: 7808x11008, 3904x5504, 1952x2752, 976x1376, 488x688

Its seems images generated with _rio_calculate_cog are sometimes bigger then the one with _gdal_calculate_cog but still tiff (1600m) -> 380m(rio) - 180(gdal) tiff (800m) -> 50m(rio) - 80(gdal)

not clear to me why.

mmaelicke commented 2 months ago

You can change either of the two functions to make them work, but please always double-Check that visualization is really working.

You need to always use openlayers as, ie QGis can show all images, even the ones openlayers fails

mmaelicke commented 2 months ago

Die unterschiedlichen Ergebnisse werden an den defaults von cogeo liegen. Mir fallen spontan die Zahl der overviews und deren Komprimierung (der ovetviews, nicht das Bild selber) ein, die unterschiedlich sein können.

JesJehle commented 2 months ago

You can change either of the two functions to make them work, but please always double-Check that visualization is really working.

You need to always use openlayers as, ie QGis can show all images, even the ones openlayers fails

it's all tested with OpenLayers in a codesanbox