Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

geoTIFF output should have band description in a different format #257

Open
clausmichele opened this issue Oct 27, 2022 · 3 comments
Open

geoTIFF output should have band description in a different format #257

clausmichele opened this issue Oct 27, 2022 · 3 comments

Comments

@clausmichele
Copy link
Member

Hi,
I'm working on the client side processing and trying to automatize the metadata extraction from netCDFs and geoTIFFs.
Currently, I'm not able to get band names from geoTIFFs generated by VITO.

You can check this issue using rasterio or gdal:

import rasterio
with rasterio.open('./results/S2_L2A_data.tiff') as ds:
    meta = ds.descriptions
print(meta)

(None, None, None)
from osgeo import gdal
rds = gdal.Open('./results/S2_L2A_data.tiff')
bands = {rds.GetRasterBand(i).GetDescription(): i for i in range(1, rds.RasterCount + 1)}
print(bands)

{'': 3}

The issue could be solved putting the band names with the field Description not in metadata (I'm not sure if it makes any difference if it's uppercase, lowercase or with the first letter capital as GDAL writes it).

A way to add the description in the "right" way is the following (S2_L2A_data_mod.tiff is a copy of S2_L2A_data.tiff):

#source https://gis.stackexchange.com/questions/416391/why-does-gdals-getdescription-function-return-an-empty-string
from osgeo import gdal
def set_raster_band_description(ras_file=None, n_band=None, s_description=None):
    """  Sets s_description to n_band of ras_file  """
    gdal_ras = gdal.Open(ras_file, gdal.GA_Update)
    gdal_ras_band = gdal_ras.GetRasterBand(n_band)
    assert not gdal_ras_band.GetDescription(), f"The raster band {n_band} already has a description"
    gdal_ras_band.SetDescription(s_description)
    assert gdal_ras_band.GetDescription(), f"The description was not set on raster band {n_band}"
    gdal_ras = None
set_raster_band_description('./results/S2_L2A_data_mod.tiff',1,'B02')
set_raster_band_description('./results/S2_L2A_data_mod.tiff',2,'B03')
set_raster_band_description('./results/S2_L2A_data_mod.tiff',3,'B04')

After this step, the band names are correctly recognized automatically:

import rasterio
with rasterio.open('./results/S2_L2A_data_mod.tiff') as ds:
    meta = ds.descriptions
print(meta)

('B02', 'B03', 'B04')

from osgeo import gdal
rds = gdal.Open('./results/S2_L2A_data_mod.tiff')
bands = {rds.GetRasterBand(i).GetDescription(): i for i in range(1, rds.RasterCount + 1)}
print(bands)

{'B02': 1, 'B03': 2, 'B04': 3}

tiff_samples.zip

gdalinfo on VITO geoTIFF 
Driver: GTiff/GeoTIFF
Files: ./results/S2_L2A_data.tiff
Size is 626, 240
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 32N",
    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["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            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["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State."],
        BBOX[0,6,84,12]],
    ID["EPSG",32632]]
Data axis to CRS axis mapping: 1,2
Origin = (660850.000000000000000,5104100.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  PROCESSING_SOFTWARE=0.6.4a1
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  660850.000, 5104100.000) ( 11d 4'47.99"E, 46d 4'17.57"N)
Lower Left  (  660850.000, 5101700.000) ( 11d 4'45.07"E, 46d 2'59.85"N)
Upper Right (  667110.000, 5104100.000) ( 11d 9'39.20"E, 46d 4'12.16"N)
Lower Right (  667110.000, 5101700.000) ( 11d 9'36.17"E, 46d 2'54.45"N)
Center      (  663980.000, 5102900.000) ( 11d 7'12.11"E, 46d 3'36.03"N)
Band 1 Block=256x256 Type=UInt16, ColorInterp=Red
  NoData Value=0
  Metadata:
    DESCRIPTION=B02
Band 2 Block=256x256 Type=UInt16, ColorInterp=Green
  NoData Value=0
  Metadata:
    DESCRIPTION=B03
Band 3 Block=256x256 Type=UInt16, ColorInterp=Blue
  NoData Value=0
  Metadata:
    DESCRIPTION=B04
gdalinfo on the modified geoTIFF 
Driver: GTiff/GeoTIFF
Files: ./results/S2_L2A_data_mod.tiff
Size is 626, 240
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 32N",
    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["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            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["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State."],
        BBOX[0,6,84,12]],
    ID["EPSG",32632]]
Data axis to CRS axis mapping: 1,2
Origin = (660850.000000000000000,5104100.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  PROCESSING_SOFTWARE=0.6.4a1
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  660850.000, 5104100.000) ( 11d 4'47.99"E, 46d 4'17.57"N)
Lower Left  (  660850.000, 5101700.000) ( 11d 4'45.07"E, 46d 2'59.85"N)
Upper Right (  667110.000, 5104100.000) ( 11d 9'39.20"E, 46d 4'12.16"N)
Lower Right (  667110.000, 5101700.000) ( 11d 9'36.17"E, 46d 2'54.45"N)
Center      (  663980.000, 5102900.000) ( 11d 7'12.11"E, 46d 3'36.03"N)
Band 1 Block=256x256 Type=UInt16, ColorInterp=Red
  Description = B02
  NoData Value=0
  Metadata:
    DESCRIPTION=B02
Band 2 Block=256x256 Type=UInt16, ColorInterp=Green
  Description = B03
  NoData Value=0
  Metadata:
    DESCRIPTION=B03
Band 3 Block=256x256 Type=UInt16, ColorInterp=Blue
  Description = B04
  NoData Value=0
  Metadata:
    DESCRIPTION=B04
@jdries
Copy link
Contributor

jdries commented Nov 7, 2022

Can you do a 'tiffinfo' on the 'correct' tiff versus the vito tiff? That will tell us how the actual tiff metadata looks like.
Note that this is a custom gdal feature, so it's actually not guaranteed at all that an openEO backend returns tiffs with this custom metadata. May be better to rely on STAC metadata instead.

@clausmichele
Copy link
Member Author

Thanks for the clarification Jeroen, I agree that even GDAL metadata fields are not the standard, but I just wanted to mention this, so that you are aware of it.

Anyway, the tiffinfo outputs are here.
original from VITO:

tiffinfo results/S2_L2A_data.tiff 
TIFFReadDirectory: Warning, Unknown field with tag 33550 (0x830e) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 33922 (0x8482) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34735 (0x87af) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 42112 (0xa480) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 42113 (0xa481) encountered.
=== TIFF directory 0 ===
TIFF Directory at offset 0x8 (8)
  Image Width: 626 Image Length: 240
  Tile Width: 256 Tile Length: 256
  Bits/Sample: 16
  Sample Format: unsigned integer
  Compression Scheme: AdobeDeflate
  Photometric Interpretation: RGB color
  Samples/Pixel: 3
  Planar Configuration: separate image planes
  Tag 33550: 10.000000,10.000000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,660850.000000,5104100.000000,0.000000
  Tag 34735: 1,1,0,5,1024,0,1,1,1025,0,1,1,2054,0,1,9102,3072,0,1,32632,3076,0,1,9001
  GDAL Metadata: <GDALMetadata>
  <Item name="PROCESSING_SOFTWARE">0.6.5a1</Item>
  <Item name="DESCRIPTION" sample="0">B02</Item>
  <Item name="DESCRIPTION" sample="1">B03</Item>
  <Item name="DESCRIPTION" sample="2">B04</Item>
</GDALMetadata>
  GDAL NoDataValue: 0
  Predictor: none 1 (0x1)

and with added GDAL description fields:

tiffinfo S2_L2A_data_mod.tiff
TIFFReadDirectory: Warning, Unknown field with tag 33550 (0x830e) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 33922 (0x8482) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34735 (0x87af) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 42112 (0xa480) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 42113 (0xa481) encountered.
=== TIFF directory 0 ===
TIFF Directory at offset 0xb5448 (742472)
  Image Width: 626 Image Length: 240
  Tile Width: 256 Tile Length: 256
  Bits/Sample: 16
  Sample Format: unsigned integer
  Compression Scheme: AdobeDeflate
  Photometric Interpretation: RGB color
  Samples/Pixel: 3
  Planar Configuration: separate image planes
  Tag 33550: 10.000000,10.000000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,660850.000000,5104100.000000,0.000000
  Tag 34735: 1,1,0,5,1024,0,1,1,1025,0,1,1,2054,0,1,9102,3072,0,1,32632,3076,0,1,9001
  GDAL Metadata: <GDALMetadata>
  <Item name="PROCESSING_SOFTWARE">0.6.4a1</Item>
  <Item name="DESCRIPTION" sample="0">B02</Item>
  <Item name="DESCRIPTION" sample="0" role="description">B02</Item>
  <Item name="DESCRIPTION" sample="1">B03</Item>
  <Item name="DESCRIPTION" sample="1" role="description">B03</Item>
  <Item name="DESCRIPTION" sample="2">B04</Item>
  <Item name="DESCRIPTION" sample="2" role="description">B04</Item>
</GDALMetadata>

  GDAL NoDataValue: 0
  Predictor: none 1 (0x1)

@jdries jdries mentioned this issue Jan 25, 2023
Open
@jdries
Copy link
Contributor

jdries commented Jan 25, 2023

Ok, it's a geotrellis issue:
locationtech/geotrellis#3496

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@jdries @clausmichele