CategoricalVariable feature.

[lupemba]

Add handling for categorical variables to CommonDataModel. This allows a variable stored as Integers to be automatically converted to a CategoricalArrays.CategoricalArray

This is possible without any updates to DiskArrays and thus makes JuliaIO/DiskArrays.jl#287 irrelevant.

The primary focus is to solve JuliaGeo/NCDatasets.jl#143 but this also requires a PR to NCDatasets.jl.

TODO

• Improve code readability
• Add documentation
• (Optional) Add write support

[lupemba]

@tiemvanderdeure and @evetion
Here is a draft PR for a CategoricalVariable. I still have to clean it up but the read part works fine for me.

[tiemvanderdeure]

Thanks - this looks much more robust than whatever I was trying in Rasters.jl.

One question is if it would be worth putting this in an extension on CategoricalArrays to avoid adding a dependency?

[tiemvanderdeure]

This would be very ineffecient - we are basically deconstructing and then reconstructing the entire categorical array. We need to pass the integer values and the mapping directly to CategoricalArray to avoid this, that will be much faster

[Alexander-Barth]

yes, it would be great if we can reuse the underling raw array and wrap it by annotating it with the labels.

This issue at CategoricalArrays.jl discusses a similar use case:
JuliaData/CategoricalArrays.jl#389

In Implementation details it is mentioned that 0 is always the missing value which can be a problem if 0 maps to something different in the NetCDF file.

[lupemba]

Yes, the CategoricalArrays.jl has some annoying limitations. This PR is not focused on performance but just on enabling the feature with a stable interface. If the interface is good then it is possible to improve the performant later on.

At the moment NCDatasets.jl throws an error for NC enum making these variables inaccessible so performance is not the first concern.

[tiemvanderdeure]

I would say the way we wrap things and convert to CategoricalArray is key to make this work properly so let's get it right now :). Especially because CategoricalArrays has quite a few limitations that we will have to hack our way around

[tiemvanderdeure]

One of the main problems is that CategoricalArrays requries that a CategoricalPool has consecutive indices starting at 1, which most files will not have. So just using the Variable itself as the underlying array in the CategoricalArray would not work. Maybe we need to do something like have a different dictionary that maps the values in the file to index values, wrap the Variable in that, and then wrap that in a CategoricalArray (where the CategoricalPool comes directly from the mapping, again with changed indices).

[Alexander-Barth]

Actually, I am wondering if improving the performance is actually possible within the bounds of the stable API of CategoricalArrays. Would we not also need to scan the whole array to handle the case when 0 is not the missing value?

But of course, the current situation (error with enums) is indeed quite bad.

[Alexander-Barth]

Hi Simon, thanks a lot for this PR! :-)
I was wondering too if it is better to have CategoricalArrays as an extension for CommonDataModel as for example Zarr does not support enums (as far as I can tell). CategoricalArrays can still be a "strong" dependency for NCDatasets so that any netcdf file works out of the box.

[tiemvanderdeure]

An alternative way to design this would be to have a concrete type that wraps an AbstractVariable - then packages could call something like CDM.categorical(var, dict) and would not have to implement their own AbstractCategoricalVariable.

I'm suspecting that these AbstractCategoricalVariables would resemble each other quite a lot so there would be redundancy there.

[Alexander-Barth]

I think there are several ways to read an array of NetCDF enums:

Using CategorialArrays: advantages is that the code is quite straightforward for arbitrary enum types.
The storage type in memory seems to be UInt32(?) while the original data can have much smaller type like Int8. It is
unclear if we can avoid copying the data array using only the public API of CategorialArrays (I would love to hear if there is a way to do that, even if it is not implemented in the PR).
Internally CategorialArrays uses 0 as missing value, but for a NetCDF enum, the value 0 can map to something else.
When writing a CategorialArrays we need to name a type and underlying storage type. The user would need to provide this information separately.

Using jula enums: The storage type can match the underlying NetCDF file and we can load the data without extra copying. We can easily query the name and storage type of the enum type when writing the data.
Creating the enum type dynamically however would involve calling eval (which is probably also necessary for mapping NetCDF combound types to julia structs). Every enum would need to be scoped inside the temporary module. The code is more complex but these tests seem to work fine. In particular, an array of enums can be directly used as a buffer for loading. Julia enums are native to julia and no additional dependency is needed. Maybe the user would be the least surprised if NetCDF enums are loaded as julia enums.

Julia array of symbols: This is a native data structure in julia. But every symbol is 8 bytes. The user would need to provide the name and storage type separately for writing.

Is there anything that I forgot or something to be corrected?

[evetion]

I've implemented flag values in HDF5 as Categorical Arrays in SpaceLiDAR here https://github.com/evetion/SpaceLiDAR.jl/blob/59c14d680ad20056c039d7f22e2f09435eb6dfca/src/H5Table/utils.jl#L636.

It seems the most pragmatic at the moment, but I like the proposal of enums backed by the original data. Makes comparison feel more logical too (to a symbol thing, not a random string). Then again, maybe its better to extend CategoricalArrays with a simple read only implementation (dict based, no compilation) that takes any backing data?

[lupemba]

One question is if it would be worth putting this in an extension on CategoricalArrays to avoid adding a dependency?

I have been a bit undecided about that. My first plan was to make an extension but then I went for the simple solution. I will turn it back into an extension.

An alternative way to design this would be to have a concrete type that wraps an AbstractVariable - then packages could call something like CDM.categorical(var, dict) and would not have to implement their own AbstractCategoricalVariable.

I'm suspecting that these AbstractCategoricalVariables would resemble each other quite a lot so there would be redundancy there.

Good idea, I will update the PR

I think there are several ways to read an array of NetCDF enums:

Using CategorialArrays: advantages is that the code is quite straightforward for arbitrary enum types. The storage type in memory seems to be UInt32(?) while the original data can have much smaller type like Int8. It is unclear if we can avoid copying the data array using only the public API of CategorialArrays (I would love to hear if there is a way to do that, even if it is not implemented in the PR). Internally CategorialArrays uses 0 as missing value, but for a NetCDF enum, the value 0 can map to something else. When writing a CategorialArrays we need to name a type and underlying storage type. The user would need to provide this information separately.

Using jula enums: The storage type can match the underlying NetCDF file and we can load the data without extra copying. We can easily query the name and storage type of the enum type when writing the data. Creating the enum type dynamically however would involve calling eval (which is probably also necessary for mapping NetCDF combound types to julia structs). Every enum would need to be scoped inside the temporary module. The code is more complex but these tests seem to work fine. In particular, an array of enums can be directly used as a buffer for loading. Julia enums are native to julia and no additional dependency is needed. Maybe the user would be the least surprised if NetCDF enums are loaded as julia enums.

Julia array of symbols: This is a native data structure in julia. But every symbol is 8 bytes. The user would need to provide the name and storage type separately for writing.

Is there anything that I forgot or something to be corrected?

CategorialArrays is a pretty common package and it seems to support the main use cases of NetCDF enums. The few examples I have seen of NetCDF enums, they are used as categorical variables. I therefore went for CategorialArrays.jl as the simplest solution to the use cases I am aware off.

Your Julia enum example looks really interesting. The reason I initially skipped the idea is, enums are Julia types. Dynamically defining types based on the content of the NetCDF adds a layer of complexity we probably don't need. I will have to try your example to see how complex it really is, but I don't see a clear advantage. Do you have a use case where you would prefer to get a "dynamically defined enum" instead of a "CategoricalArray"?

The Julia enum solution might turn out to be elegant, but the CategoricalArray solution is so simple that I can both implement and use it 😅

[Alexander-Barth]

Do you have a use case where you would prefer to get a "dynamically defined enum" instead of a "CategoricalArray"?

One use case I can think of is when you copy a variable from one NetCDF file to another. An array of Julia enums contains all the information (type name, underlying storage and class to integer mapping, and order of enums) to create the same data. The mapping integers <-> String of CategoricalArray is not the same as the enum's defined on the NetCDF file. For data that you create for your self, this might not be important, but for data that you distribute publicly it might be preferable to avoid unnecessary variations. It is possible to have stable integer mapping with CategoricalArray but additional work would be necessary (like copying the data).

Another aspect could be loading data in-place (NCDatasets.load!) to avoid allocations of large arrays. For enums, it approach is quite straightforward as the array of enums can be passed to the C-layer. For CategoricalArray we would (probably?) need to use the internal API.

I made some tests with CategoricalArrays.jl and apparently Symbols are not supported.
Unfortunately, comparing strings is not very fast:

julia> @btime count(==("Clear"),data5) # CategoricalVector
  520.717 μs (3 allocations: 48 bytes)
8375
julia> @btime count(==(mod.Clear),data4) # Vector of enums
  8.912 μs (3 allocations: 48 bytes)
8375

julia> typeof(data5)
CategoricalVector{String, Int8, String, CategoricalValue{String, Int8}, Union{}} (alias for CategoricalArray{String, 1, Int8, String, CategoricalValue{String, Int8}, Union{
}})

julia> typeof(data4)
Vector{cloud_class_t} (alias for Array{Main.mod_cloud_class_t.cloud_class_t, 1})

julia> length(data4), length(data5)
(100000, 100000)

I made some benchmarks, to read these 100 000 elements, and these are the run time that I got:

# loading enum
  17.551 μs (5 allocations: 97.79 KiB)                                                                                                                                       
# loading and converting to CategoricalVector
  3.154 ms (50 allocations: 978.98 KiB) 

(see @tiemvanderdeure comments) . However, I never used CategoricalArrays before so let me know if I did something stupid 😅. I know that we can always optimize later, but I don't have the impression that manipulating directly the refs field of a CategoricalArray is public API.

But maybe there are other use-cases that I am not aware of where CategoricalArrays is better suited than vector of enums 😀.

[lupemba]

@Alexander-Barth
Maybe I focused to much on the "read" part and forgot the "write" part. When the feature is extend to have write support enums probably fell more natural. I still have a lot of open question in my head.

• Does the dynamic definition of enum trigger a lot of compilation. Julia will need to compile new functions for the new enum types.
• If we read multiple files with the same enum will they have the same type or do we define one per file/varible.
• Can the user define the enum up front and then use that type to read the variable, for example in scripts where you both write and read enum varibles.
• Does AbstractDiskArray{<:Enum} work as expected or will we have to make updates to DiskArray.jl

It will probably be possible to handle a lot of the higher level dispatch with ::AbstractVariable{<:Enum}

[Alexander-Barth]

@lupemba these are excellent questions and actually while working compound types I came across these as well 😀.

• In my test compiling the example enum took 50ms in a fresh julia session. It is indeed much larger than the IO for 100_000 values, but you pay this only once. In the case of compound types, the struct are cached at the level of data sets. We might use the same of enums.

• "same type or do we define one per file/variable": For compound types, I currently implemented that the reconstructed struct should be the same of the entire file as they are globally defined similar to dimensions. But they are different (per default) for different files. But the user can easily query the type to force them to be the same.

• For compound types, the user can force NCDatasets to reuse MyStruct for the netCDF custom type nc_compound_t with the usertype! method (defined in the branch https://github.com/JuliaGeo/NCDatasets.jl/tree/ab-compound-types):

usertype!(ds,"nc_compound_t",MyStruct)

We can also query the type to be used without reading any data with usertype(ds,"nc_compound_t"). It will return the reconstructed type unless the user already set the type before. That mechanism can also be used for enums if the user has already the enum defined.

Whether all this works with DiskArrays of course remains to be seen. But al least with an array of custom struct it did work without any problems.

Another related point, that I came across, is that NetCDF enums can also be fields of user compound types:

netcdf weather_data {
types:
  int enum cloud_t {Clear = 0, Cumulonimbus = 1, Stratus = 2, Overcast = 3};

  compound obs_t {
    float temperature;
    int humidity;
    cloud_t sky_condition;  // Nested Enum
  };

dimensions:
  station = 3;

variables:
  obs_t weather_reports(station);
...
}

With the corresponding definition of the struct and enum, nc_put_var and nc_get_var! can already handle such arrays in my tests. If we would use CategoricalArrays, these fields would need to be converted "manually" to CategoricalValue (and needing to copy the data of the whole loaded array as far as I can see).

[lupemba]

@Alexander-Barth
Okay I am convinced now.
I am still a bit uncertain, what logic is general (to be placed in CommonDataModel) and what logic is NetCDF specific?

@tiemvanderdeure, @evetion
Will arrays of Julia Enums also fit your use cases?

[tiemvanderdeure]

The only thing I would say in favour of doing this with CategoricalArrays is that other packages use it as well. So if a Raster with land cover classes is an input for a machine learning model with MLJ, it will "just work" if it is read as a CategoricalArray, whereas with an enum I need to do some kind of conversion.

Maybe we can solve this by making this conversion easy, but it's still an extra step.

I did a re-write of this branch in lupemba#3, where reading is much faster and CategoricalVariable is a concrete type that contains any AbstractVariable and a mapping.

[Alexander-Barth]

@lupemba Actually, I only needed changes in NCDatasets to support structs (netcdf compound types). The container array is unchanged just the element type is special. I can easily make an alternative PR with enums to test this.

@tiemvanderdeure
Yes, I think it would be nice if CategoricalArray could take an arrays of enums. These array types are so similar, I was a bit surprised that this does not work already.

As far as I can tell, the .pool field is private for CategoricalArray (as it is mention in Implementation-details) so I would avoid to access it directly in CommonDataModel.jl, but of course CategoricalArray, when receiving an arrays of enums, can manipulate it directly.