NetCDF to CoG workflow

[CapeHobbit]

Hi cogeotiff team,

I apologise in advance for this very general query. I am looking for pythonic ways to translate level-2 satellite data from a series of 2D irregularly gridded netCDF files to CoG. I very much see rio-cogeo as the last point in this chain, but I am wondering if you have recommendations about how to reach to a good "ingestion point" for it. So far, reprojecting from the original grid with rasterio/rio-xarray is causing me huge headaches (and yet they seem to recommended over osgeo gdal), and so I wondered if you have encountered, or can recommend, a similar workflow?

Things I want to avoid;

• using large external packages (e.g. SNAP)
• using GPL3 software (e.g. Satpy)

I realise that my problem is not the one you are trying to solve, but any advice would be appreciated.

[vincentsarago]

@CapeHobbit

You have multiple options here. You could use GDAL python bindings to do everything or rasterio+rio-cogeo

You could in theory do the reprojection/conversion to COG in one GO (without needed to write intermediate files) but it will depends on how your original NetCDF is structured.

import rasterio
from rasterio.vrt import WarpedVRT
from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles

out_profile = cog_profiles.get("deflate")

with rasterio.open("my.nc") as src:
    with WarpedVRT(src, crs="epsg:3857") as vrt:
        cog_translate(vrt, "out.tif", out_profile, allow_intermediate_compression=True)

[CapeHobbit]

@vincentsarago

Thanks for the quick reply and suggested code. I would love to do it this way, but unfortunately the data I am working with (Sentinel-3 OLCI) separates the bands (one per netCDF file) and lat/lon (a separate netcdf file). The projection information (such that it is) is held in an associated XML file. When I run your suggested code, I get the following error;

Unable to compute a transformation between pixel/line and georeferenced coordinates for ./test_data/S3B_OL_2_WFR____20230507T094601_20230507T094901_20230508T205812_0179_079_136_1980_MAR_O_NT_003.SEN3/geo_coordinates.nc. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

I tried with both epsg:3857 and epsg:4326, with the same result - which is not surprising as there is no defined projection, just lat/lon WGS84 coords! I independently made some small progress using pyresample to reduce the coordinates to 1 dimension. I can write this to a tiff with the following spatial_ref, but when I try to convert it to geotiff, I get the following error;

Cannot invert geotransform

I think this means that I have not defined the transform properly in the first place. My failing test tiff has the following spatial_ref;

<xarray.DataArray 'spatial_ref' ()>
array(0)
Coordinates:
    spatial_ref  int64 0
Attributes:
    crs_wkt:                      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...
    semi_major_axis:              6378137.0
    semi_minor_axis:              6356752.314245179
    inverse_flattening:           298.257223563
    reference_ellipsoid_name:     WGS 84
    longitude_of_prime_meridian:  0.0
    prime_meridian_name:          Greenwich
    geographic_crs_name:          WGS 84
    horizontal_datum_name:        World Geodetic System 1984
    grid_mapping_name:            latitude_longitude
    spatial_ref:                  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...

Projections/CRS etc are not my strong point, so thanks in advance if you can provide any more help here. Below is a link to the netCDF files and XML I referenced above.

Sentinel3_files.zip

[CapeHobbit]

As a quick update; I believe I have managed to make some progress. If I use pyresample to convert the 2D lat/lon grid to 1D coords projected in epsg:3857, then I can parse this data into a new xarray dataset, save as a tiff and then convert to a cogeotiff. Barring some small issues with transparency, the eventual CoG renders properly in QGIS. I am not sure why epsg:4326 does not work, but I will investigate. My solution is a little longwinded, using pyresample, pyproj, rioxarray and cogeo but no explicit rasterio. I can share it here, if it is of interest.

There are about 100 stack overflow posts on this, none with satisfactory answers (at least for me!) so perhaps refining this method would be of use?

[vincentsarago]

Here is a solution I though about but it doesn't seem to give a correct result 🤔

import numpy
import netCDF4 as nc

from rasterio.transform import from_bounds

ds = nc.Dataset("chl_nn.nc")
coord = nc.Dataset("geo_coordinates.nc")

lat = coord["latitude"][:]
lon = coord["longitude"][:]

xmin = numpy.amin(lon)
xmax = numpy.amax(lon)
ymin = numpy.amin(lat)
ymax = numpy.amax(lat)

print(xmin, ymin, xmax, ymax)

height, width = ds["CHL_NN"].shape
transform = from_bounds(xmin, ymin, xmax, ymax, width, height)
print(transform)


import rasterio
from rio_tiler.utils import _array_gdal_name
from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles

data = ds["CHL_NN"][:, :]
datasetname = _array_gdal_name(data)

with rasterio.open(datasetname, "r+") as src:
    src.transform = transform
    src.crs = "EPSG:4326"
    prof = cog_profiles.get("deflate")
    cog_translate(src, "cog.tif", prof)

[CapeHobbit]

Ah, I know what the problem may be! I found that the latitude array in this granule is somehow "corrupted". There are a couple of scan line where the latitude is missing or zero, which throws the whole things off! Simple(-ish) solution is to choose your max/min from the image corners (not trustable in all cases, as some swaths have bends that are wider in extent than the corners, but i'll deal with that later!). Your solution is about 5 times shorter than mine! :)