Ingesting simple gridded NetCDF data

[petwa]

I have managed to set  up a Syntool installation and I would now like to add some of our data. I have looked at the Developer Guide and the reader_skeleton.py but I would find it very useful to see an example that closely aligned to our data. 

We already have our data mapped  in NetCDF files so the conversion to a Syntool Geotiff should not require much work. The reader_skeleton.py does not seem to be focused on this usecase; is there a more suitable example converter I could use as a basis? Is there any other source of documentation for this other than the Developer Guide?

 

Cheers, 

Pete 

[Sylvain Herlédan]

Hi Pete,

I think the developer guide and other manuals available on https://ftp.odl.bzh/odl/syntool/ are the only documentation we published about Syntool.

In most cases there are three product-specific steps that you have to look out for when writing a new Syntool reader:

• choosing the mask to apply on the data variables (requires data expertise to make sensible choices)
  
  
• defining geolocation so that the resulting GeoTIFF can be warped by GDAL with sufficient precision during the ingestion step. For L3 and L4 products it is usually a matter of defining the geotransform (a 2x3 transformation matrix, flattened in a 6-items list in Syntool readers) and the coordinates system in which the grid is regular (lat/lon or y/x with a well-known projection) whereas for L1 and L2 data you have to use GCPs ("Ground Control Points" that associate the indices of a subset of cells in the data matrix with their geographical coordinates) because the data grid uses the instrument geometry and not a well-known projection.
  
  
• choosing the min and max values to preserve data dynamics as much as possible when packing floating point geophysical  values into 8 bits GeoTIFF bands (i.e. only 255 possible values because value '255' is reserved for transparency) .
  

The syntool_converter package contains several readers for data defined on a regular lat/lon grid and stored as NetCDF files: you can have a look at the syntool_converter/smos/smos_l4_sss_ifr.py file in the converter source code for a basic example. If it does not match your use case, could you provide more information about the data you want to add so I can point you to a more relevant example?

Cheers,

Sylvain

[petwa]

Thanks Sylvain.
I ended up using the bathymetry_gebco.py as the very simple example and hacking out large chunks of the skeleton. On visual inspection that seemed to work OK.
My next task is to get it to use our colormap. It looks like the maps are hardcoded. Is there any way to use custom palettes without hacking the syntool code?

[petwa]

Hi Sylvain,

Thanks for the syntool_converter/smos/smos_l4_sss_ifr.py, it matched pretty well and confirmed I was heading in the right direction.

Cheers,
Pete

[Sylvain Herlédan]

Hi Pete,

In most of our readers we use a helper method (see format_colortable and load_colormap defined in syntool_converter/utils/syntoolformat.py) that supports matplotlib and a few other colormaps to build the object that we assign to the colortable field of the band.

Note that using the helper method is not mandatory:  as long as a gdal.Colortable instance (wherever it comes from) matching your colormap is assigned to the colortable field, the resulting GeoTIFF will use a palette based on your colormap to render unsigned byte values as RGB.

If your colormap is really not supported (i.e. not a matplotlib colormap nor one of the other colormaps supported by Syntool), I think you should make a copy of the format_colortable method in your reader and implement a custom load_colormap method which builds and returns a matplotlib.colors.LinearSegmentedColormap for your colormap.

Cheers,

Sylvain

[petwa]

Hi Sylvain,
I was looking at load_colormap as it was used in my skeleton converter. I wanted to use the viridis colormap but it looked like the function only allowed selection from a hardcoded list and raised an exception otherwise. Is there a newer version that includes the later matplotlib colormaps?

Cheers,
Pete

[Sylvain Herlédan]

Hi Pete,

My bad, in the version available on our website the converter only supports a few matplotlib colormaps, whereas we support all of them in the development version.

If you feel comfortable with using the undocumented, unsupported dev version you can find it in the attachments, otherwise you can use the workaround I suggested in my previous message by including this code in your reader:

Code:

import gdal
import numpy as np
import matplotlib.cm
import syntool_converter.utils.syntoolformat as stfmt

def load_colormap(name):
    """
    """
    if name.startswith('matplotlib_'):
        _, cmap_name = name.split('_', 1)
        colormap = getattr(matplotlib.cm, cmap_name)
        nodata_color = (255, 255, 255)
        return colormap, nodata_color
    return stfmt.load_colormap(name)


def format_colortable(name, vmin=0., vmax=1., vmin_pal=0., vmax_pal=1.,
                      index_min=0, index_max=254, index_nodata=255):
    """
    """
    colormap, nodata_color = load_colormap(name)
    norm_min = (vmin - vmin_pal) / float((vmax_pal - vmin_pal))
    norm_max = (vmax - vmin_pal) / float((vmax_pal - vmin_pal))
    ncols = int(index_max) - int(index_min) + 1
    colors = colormap(np.linspace(norm_min, norm_max, num=ncols))
    colors = np.round(colors*255)
    entries = [(int(c[0]), int(c[1]), int(c[2])) for c in colors]
    colortable = gdal.ColorTable()
    for index in range(int(index_min), int(index_max)+1):
        colortable.SetColorEntry(index, entries[index-int(index_min)])
    if index_nodata != None:
        colortable.SetColorEntry(index_nodata, nodata_color)
    return colortable

And then pass 'matplotlib_viridis' as first argument to format_colortable.

Cheers,

Sylvain

[petwa]

Hi Sylvain,
Thanks for the help; much appreciated.
I had taken a local copy of the format_colortable and am using that. Is load_colormap used anywhere outside of format_colortable as I was lazy and set the colormap explicitly in there as we will only use the matplot colormaps?

Cheers,
Pete

[Sylvain Herlédan]

load_colormap is used in two places:

• format_colortable which is used by many readers
  
• a tool which generates a picture of the colormap that can be used as a legend in the portal. Note that it is not mandatory to generate this picture, that you can generate this image independently (with a drawing tool for example) and that we don't use these pictures anymore in the dev version of Syntool (the portal reads min, max and colormap from the config to build the legend)
  

So having a local version of format_colortable with the hardcoded colormap should not be an issue.

[petwa]

Is there a public facing Git repo that contains the latest code? I noticed there are a lot of features on the OceanDataLab sites that I can't find in the released version.

[Sylvain Herlédan]

The latest code is only available on our internal server for now because we are still fixing minor bugs and the documentation is not ready yet. We still need a few weeks to tidy things up so I think there should be a new release during the summer. We can send the source code of the dev version to users who request it but we will only provide support for the official release.