Convert NetCDF to COG modifies the Corner Coordinates

[TheRed86]

Hi,

I am using rio-cogeo to translate a netcdf file to COG.

I am working with global maps of ERA5 products, which have the following bounds: -180.0, 90.0, 179.75, -90.0

When I convert the data to COG, the Corner Coordinates I see inspecting the file (rio cogeo info file_name) are the following:

Driver: GTiff
File: /Volumes/aliastore/ERA5_exploration/t2m_new_crs/out.tif
COG: True
Compression: DEFLATE
ColorSpace: None

Profile
    Width:            1440
    Height:           721
    Bands:            53
    Tiled:            True
    Dtype:            float64
    NoData:           nan
    Alpha Band:       False
    Internal Mask:    False
    Interleave:       PIXEL
    ColorMap:         False
    ColorInterp:      ('gray', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined', 'undefined')
    Scales:           (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
    Offsets:          (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

Geo
    Crs:              EPSG:4326
    Origin:           (-180.125, 90.125)
    Resolution:       (0.25, -0.25)
    BoundingBox:      (-180.125, -90.125, 179.875, 90.125)
    MinZoom:          0
    MaxZoom:          2

Why is the origin in (-180.125, 90.125) and bounding box (-180.125, -90.125, 179.875, 90.125)?
How can I set the bounding box to be the original ones? I.e., (-180.0, 90.0, 179.75, -90.0)?

Any help would be appreciated!
Regards,
Monica

[vincentsarago]

@TheRed86 can you run gdalinfo or rio cogeo info on your input file and post the result here ?

[TheRed86]

Here you are:

 gdalinfo ../test_bounds/t2m_19400107.nc
Driver: netCDF/Network Common Data Format
Files: ../test_bounds/t2m_19400107.nc
Size is 1440, 721
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    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]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.125000000000000,90.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  latitude#axis=Y
  latitude#long_name=latitude
  latitude#standard_name=latitude
  latitude#units=degrees_north
  latitude#_FillValue=nan
  longitude#axis=X
  longitude#long_name=longitude
  longitude#standard_name=longitude
  longitude#units=degrees_east
  longitude#_FillValue=nan
  spatial_ref#crs_wkt=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  spatial_ref#geographic_crs_name=WGS 84
  spatial_ref#grid_mapping_name=latitude_longitude
  spatial_ref#horizontal_datum_name=World Geodetic System 1984
  spatial_ref#inverse_flattening=298.257223563
  spatial_ref#longitude_of_prime_meridian=0
  spatial_ref#prime_meridian_name=Greenwich
  spatial_ref#reference_ellipsoid_name=WGS 84
  spatial_ref#semi_major_axis=6378137
  spatial_ref#semi_minor_axis=6356752.314245179
  spatial_ref#spatial_ref=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  t2m#coordinates=time
  t2m#grid_mapping=spatial_ref
  t2m#long_name=2 metre temperature
  t2m#units=K
  t2m#_FillValue=nan
Corner Coordinates:
Upper Left  (-180.1250000,  90.1250000) (180d 7'30.00"W, 90d 7'30.00"N)
Lower Left  (-180.1250000, -90.1250000) (180d 7'30.00"W, 90d 7'30.00"S)
Upper Right ( 179.8750000,  90.1250000) (179d52'30.00"E, 90d 7'30.00"N)
Lower Right ( 179.8750000, -90.1250000) (179d52'30.00"E, 90d 7'30.00"S)
Center      (  -0.1250000,   0.0000000) (  0d 7'30.00"W,  0d 0' 0.01"N)
Band 1 Block=1440x1 Type=Float32, ColorInterp=Undefined
  NoData Value=nan
  Unit Type: K
  Metadata:
    coordinates=time
    grid_mapping=spatial_ref
    long_name=2 metre temperature
    NETCDF_VARNAME=t2m
    units=K
    _FillValue=nan

I see that the gdalinfo on the NetCDF file reports the same bounds but when I open the file (with xarray or Panoply) I see the original bounds in there:

file = '/Volumes/aliastore/ERA5_exploration/test_bounds/t2m_19400107.nc'

ds = xr.open_dataset(file)
ulx =  ds['longitude'].values.min() #georeferenced longitude minimum
uly =  ds['latitude'].values.max() #georeferenced latitude maximum
lrx =  ds['longitude'].values.max() #georeferenced longitude maximum
lry =  ds['latitude'].values.min() #georeferenced latitude minimum

print(ulx, uly, lrx, lry)

--> -180.0 90.0 179.75 -90.0

Panoply:

[vincentsarago]

👀

 gdalinfo ../test_bounds/t2m_19400107.nc
...
Origin = (-180.125000000000000,90.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)

The way rio-xarray / xarray and gdal calculate the bounds is slightly different corteva/rioxarray#645

[TheRed86]

Thank you very much for explaining it to me!

I understand that is not a problem with the conversion in COG but instead with the original data :(

Is there a way to fix this problem? After processing, I load the data on Geoserver, and in fact, I see the coordinates to be calculated from the LeftUpper corner...

What about using some sort of "re-bounding"? such as

gdal_translate -of COG -a_srs EPSG:4326 -a_ullr ulx, uly, lrx, lry NETCDF:"data.nc":var out.tiff

Is there a way to pass the option -a_ullr to your wrapper library?

[TheRed86]

Dear @vincentsarago,

I apologise if my question was stupid.
Unfortunately, I am pretty new to COG...

I need to learn as much as I can to work with COG in Cloud.

Assume I retrieve an hourly dataset from a provider (for example, an ERA5 dataset from ECMWF).
From this dataset, I need to compute the weekly average and serve on a web application each weekly value as a COG file, creating an ImageMosaic in Geoserver. The COG files must be stored on an Object Storage S3.

I see some helpful functionality on rio-cogeo that could help me to achieve this.

For now, I have used rio-cogeo in this way to create the COG files with the weekly maps:

def calculate_weekly_means_and_save_cog(file, var, output_dir, output_file):
    print(" * Calculating weekly means ... ")
    ds = xr.open_dataset(file, engine='netcdf4')
    ds = ds.sortby('time', ascending=True)

    # Resample the data to weekly frequency and calculate the mean
    weekly_means = ds[var].resample(time='1W', closed='left', label='right').mean(dim='time')
    weekly_means.to_netcdf(output_file)

    # Create the output directory if it doesn't exist
    os.makedirs(output_dir, exist_ok=True)

    # Loop through the weeks and save the weekly means to separate COG files
    for i, time in enumerate(weekly_means.time.values):
        week_start = pd.to_datetime(time).strftime('%Y%m%d')
        output_name = os.path.join(output_dir, f't2m_{week_start}.nc')
        weekly_mean = weekly_means[i]
        cog_file = output_name[:-2] + "tif"
        wm_4326 = weekly_mean.rio.write_crs("EPSG:4326")
        wm_4326.to_netcdf(output_name)
        command = [
            "rio", "cogeo", "create", output_name, cog_file
        ]
        subprocess.run(command)
        os.remove(output_name)

    ds.close()
    weekly_means.close()
    return "Weekly means computed and saved to separate files!"

Is there a better way, using the Python API (https://cogeotiff.github.io/rio-cogeo/API/), to accomplish the COGs generation?
What about "Create output COG in Memory"? I see that it is designed to upload COG in S3.

Best regards!

[vincentsarago]

I apologise if my question was stupid.

ah ah no, there is no stupid question, I just didn't had the time to answer.

You do not need to crop your file, because I guess you want to keep your output COG as your input dataset.
The dataset bounds is not wrong even if they are a bit weird!

as explained in corteva/rioxarray#645 your input file (netcdf) is a grid and the coordinates should be the centroid while it seems they represent the upper left value, which is why when converting to GeoTIFF the cordinates a slightly different (because bounds are upper left pixel coordinates).

Is there a better way, using the Python API (https://cogeotiff.github.io/rio-cogeo/API/), to accomplish the COGs generation?
What about "Create output COG in Memory"? I see that it is designed to upload COG in S3.

Yeah for your usecase, using the direct python API will be much simpler.

[TheRed86]

Thank you very much for your patience!

You are right. I initially have global data and want the output to remain a global map. However, thanks for clarifying.

Yeah for your usecase, using the direct python API will be much simpler.

I tried, but I got an error :(
The idea is not to generate an intermediate NetCDF file containing the single timestep that I have to convert to COG but work with the xarray.dataset directly.
I tried to modify my function in the following way:

def calculate_weekly_means_and_save_cog(file, var, output_dir, output_file):
    print(" * Calculating weekly means ... ")
    with xr.open_dataset(file, engine='netcdf4') as ds:
        ds = ds.sortby('time', ascending=True)

        # Resample the data to weekly frequency and calculate the mean
        weekly_means = ds[var].resample(time='1W', closed='left', label='right').mean(dim='time')

        # Create the output directory if it doesn't exist
        os.makedirs(output_dir, exist_ok=True)

        # Loop through the weeks and save the weekly means to separate COG files
        for i, time in enumerate(weekly_means.time.values):
            week_start = pd.to_datetime(time).strftime('%Y%m%d')
            output_name = os.path.join(output_dir, f't2m_{week_start}.tif')
            weekly_mean = weekly_means[i]

            # Prepare additional arguments for COG creation
            dst_kwargs = {
                "compress": "deflate",  # Compression type
            }
            # Convert to COG directly without intermediate NetCDF file
            cog_translate(weekly_mean, output_name, in_memory=False, dst_kwargs=dst_kwargs)

    return "Weekly means computed and saved to separate COG files!"

but I get the following error:

File "/Users/thered86/miniconda3/envs/era5_env/lib/python3.12/site-packages/rio_cogeo/cogeo.py", line 181, in cog_translate
    src_dst = ctx.enter_context(rasterio.open(source))
                                ^^^^^^^^^^^^^^^^^^^^^
  File "/Users/thered86/miniconda3/envs/era5_env/lib/python3.12/site-packages/rasterio/env.py", line 451, in wrapper
    return f(*args, **kwds)
           ^^^^^^^^^^^^^^^^
  File "/Users/thered86/miniconda3/envs/era5_env/lib/python3.12/site-packages/rasterio/__init__.py", line 216, in open
    raise TypeError("invalid path or file: {0!r}".format(fp))
TypeError: invalid path or file: <xarray.DataArray 't2m' (latitude: 721, longitude: 1440)>

I apologise, but I don't understand from the documentation how to use the API without saving the intermediate NetCDF file with the single step.

If I do:

 for i, time in enumerate(weekly_means.time.values):
            week_start = pd.to_datetime(time).strftime('%Y%m%d')
            output_name_nc = os.path.join(output_dir, f't2m_{week_start}.nc')
            output_name_tif = os.path.join(output_dir, f't2m_{week_start}.tif')
            weekly_mean = weekly_means[i]
            weekly_mean.to_netcdf(output_name_nc)

            out_profile = cog_profiles.get("deflate")
            # Convert to COG directly without intermediate NetCDF file
            with rasterio.open(output_name_nc) as src:
                with WarpedVRT(src, crs="epsg:3857") as vrt:
                    cog_translate(vrt, output_name_tif, out_profile, allow_intermediate_compression=True)
            os.remove(output_name_nc)

Then it work, but I in fact saved the intermediate file :(