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.

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

Store band specific attributes in coordinates #81

Closed
snowman2 opened this issue Jan 21, 2020 · 3 comments · Fixed by #600
Closed

Store band specific attributes in coordinates #81

snowman2 opened this issue Jan 21, 2020 · 3 comments · Fixed by #600
Labels
proposal Idea for a new feature.
Milestone

Comments

@snowman2
Copy link
Member

snowman2 commented Jan 21, 2020

Currently there are band-specific attributes stored in the attrs. However, this can be potentially lossy when performing xarray operations. So, I am thinking about storing them in the coordinates as arrays along the band dimension. This will also be useful when selecting a specific band as you will only have the subset of data associated with the band.

Proposed coordinates:

  • nodatavals
  • scales
  • offsets
  • descriptions (maybe)
  • units (maybe)

Additionally. I am thinking about adding a profile and or metadata coordinate to store the profile information from the raster not stored in other attributes to persist when writing. For example - compression, tiling, etc...

This should make things simpler and more reliable for #76.

cc: @alfredoahds @justingruca

@snowman2
Copy link
Member Author

Not much point in adding nodatavals: rasterio/rasterio#1859

@clausmichele
Copy link

Is this related to the issue I'm seeing? Loading a geoTIFF with multiple bands, with their names stored in each band metadata, only the first one is loaded in the attrs:

rioxarray output:

print(rioxarray.open_rasterio('results/S2_L2A_data.tiff'))

<xarray.DataArray (band: 3, y: 240, x: 626)>
[450720 values with dtype=uint16]
Coordinates:
  * band         (band) int64 1 2 3
  * x            (x) float64 6.609e+05 6.609e+05 ... 6.671e+05 6.671e+05
  * y            (y) float64 5.104e+06 5.104e+06 ... 5.102e+06 5.102e+06
    spatial_ref  int64 0
Attributes:
    DESCRIPTION:   B02
    _FillValue:    0.0
    scale_factor:  1.0
    add_offset:    0.0

gdalinfo output:

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

I attach the geoTIFF I'm using for further tests.

S2_L2A_data.zip

@snowman2
Copy link
Member Author

Yes, it is also related to #296

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
proposal Idea for a new feature.
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants