Iris not able to read certain variables on a file

[Bas Crezee]

Hi all, 

I am using Iris to read landcover data recently released on the Climate Data Store. In case someone wants to reproduce the problem, I provide instructions for getting the data here. Below I copied the download script using the (Python) CDS api which I recommend, alternatively one can use the web form as well.

import cdsapi
import os


# NB: even though we request tar.gz, the downloaded file (as of 30 Oct 2019)
# is netcdf and should be renamed to have extension .nc




# Set parameters
year = 2015
outdir = '/XXX/XXXX/landcover/'


c = cdsapi.Client()


filepath = os.path.join(outdir, 'landcover_{0}.nc'.format(year))


c.retrieve(
    'satellite-land-cover',
    {
        'variable':'all',
        'format':'tgz',
        'version':'v2.0.7cds',
        'year': str(year)
    },
    filepath)

I checked the file with the cfchecker and the format is largely OK. For completion I put the output here:

CHECKING NetCDF FILE: landcover_2015.nc
=====================
Using CF Checker Version 3.2.0rc1
Checking against CF Version CF-1.6
Using Standard Name Table Version 69 (2019-10-15T14:39:18Z)
Using Area Type Table Version 9 (07 August 2018)
Using Standardized Region Name Table Version 4 (18 December 2018)




------------------
Checking variable: lccs_class
------------------


------------------
Checking variable: processed_flag
------------------


------------------
Checking variable: current_pixel_state
------------------


------------------
Checking variable: observation_count
------------------
INFO: (3.1): No units attribute set.  Please consider adding a units attribute for completeness.


------------------
Checking variable: change_count
------------------
INFO: (3.1): No units attribute set.  Please consider adding a units attribute for completeness.


------------------
Checking variable: lat
------------------


------------------
Checking variable: lon
------------------


------------------
Checking variable: crs
------------------
WARN: (3): No standard_name or long_name attribute specified


------------------
Checking variable: lat_bounds
------------------


------------------
Checking variable: lon_bounds
------------------


------------------
Checking variable: time_bounds
------------------


------------------
Checking variable: time
------------------


ERRORS detected: 0
WARNINGS given: 1
INFORMATION messages: 2

Now we come to my actual problem, it only reads in two variables properly, the others are simply neglected:

import iris

infile = '/net/exo/landclim/PROJECTS/C3S/datadir/play/landcover/landcover_2015.nc'

iris.load(infile)

cubelist = iris.load(infile)

print(cubelist)

Output:

0: crs / (1)                           (scalar cube)

1: land_cover_lccs / (1)               (time: 1; latitude: 64800; longitude: 129600)

For completion, also an ncdump of the header: 

netcdf landcover_2015 {

dimensions:

lat = 64800 ;

lon = 129600 ;

time = 1 ;

bounds = 2 ;

variables:

ubyte lccs_class(time, lat, lon) ;

lccs_class:standard_name = "land_cover_lccs" ;

lccs_class:flag_colors = "#ffff64 #ffff64 #ffff00 #aaf0f0 #dcf064 #c8c864 #006400 #00a000 #00a000 #aac800 #003c00 #003c00 #005000 #285000 #285000 #286400 #788200 #8ca000 #be9600 #966400 #966400 #966400 #ffb432 #ffdcd2 #ffebaf #ffc864 #ffd278 #ffebaf #00785a #009678 #00dc82 #c31400 #fff5d7 #dcdcdc #fff5d7 #0046c8 #ffffff" ;

lccs_class:long_name = "Land cover class defined in LCCS" ;

lccs_class:valid_min = 1 ;

lccs_class:valid_max = 220 ;

lccs_class:ancillary_variables = "processed_flag current_pixel_state observation_count change_count" ;

lccs_class:flag_meanings = "no_data cropland_rainfed cropland_rainfed_herbaceous_cover cropland_rainfed_tree_or_shrub_cover cropland_irrigated mosaic_cropland mosaic_natural_vegetation tree_broadleaved_evergreen_closed_to_open tree_broadleaved_deciduous_closed_to_open tree_broadleaved_deciduous_closed tree_broadleaved_deciduous_open tree_needleleaved_evergreen_closed_to_open tree_needleleaved_evergreen_closed tree_needleleaved_evergreen_open tree_needleleaved_deciduous_closed_to_open tree_needleleaved_deciduous_closed tree_needleleaved_deciduous_open tree_mixed mosaic_tree_and_shrub mosaic_herbaceous shrubland shrubland_evergreen shrubland_deciduous grassland lichens_and_mosses sparse_vegetation sparse_tree sparse_shrub sparse_herbaceous tree_cover_flooded_fresh_or_brakish_water tree_cover_flooded_saline_water shrub_or_herbaceous_cover_flooded urban bare_areas bare_areas_consolidated bare_areas_unconsolidated water snow_and_ice" ;

lccs_class:flag_values = 0UB, 10UB, 11UB, 12UB, 20UB, 30UB, 40UB, 50UB, 60UB, 61UB, 62UB, 70UB, 71UB, 72UB, 80UB, 81UB, 82UB, 90UB, 100UB, 110UB, 120UB, 121UB, 122UB, 130UB, 140UB, 150UB, 151UB, 152UB, 153UB, 160UB, 170UB, 180UB, 190UB, 200UB, 201UB, 202UB, 210UB, 220UB ;

byte processed_flag(time, lat, lon) ;

processed_flag:long_name = "LC map processed area flag" ;

processed_flag:standard_name = "land_cover_lccs status_flag" ;

processed_flag:valid_min = 0 ;

processed_flag:valid_max = 1 ;

processed_flag:_FillValue = -1b ;

processed_flag:flag_meanings = "not_processed processed" ;

processed_flag:flag_values = 0b, 1b ;

byte current_pixel_state(time, lat, lon) ;

current_pixel_state:long_name = "LC pixel type mask" ;

current_pixel_state:standard_name = "land_cover_lccs status_flag" ;

current_pixel_state:valid_min = 0 ;

current_pixel_state:valid_max = 5 ;

current_pixel_state:_FillValue = -1b ;

current_pixel_state:flag_meanings = "invalid clear_land clear_water clear_snow_ice cloud cloud_shadow" ;

current_pixel_state:flag_values = 0b, 1b, 2b, 3b, 4b, 5b ;

ushort observation_count(time, lat, lon) ;

observation_count:long_name = "number of valid observations" ;

observation_count:standard_name = "land_cover_lccs number_of_observations" ;

observation_count:valid_min = 0 ;

observation_count:valid_max = 32767 ;

ubyte change_count(time, lat, lon) ;

change_count:long_name = "number of class changes" ;

change_count:valid_min = 0 ;

change_count:valid_max = 100 ;

double lat(lat) ;

lat:units = "degrees_north" ;

lat:long_name = "latitude" ;

lat:standard_name = "latitude" ;

lat:valid_min = -90. ;

lat:valid_max = 90. ;

lat:bounds = "lat_bounds" ;

lat:axis = "Y" ;

double lon(lon) ;

lon:units = "degrees_east" ;

lon:long_name = "longitude" ;

lon:standard_name = "longitude" ;

lon:valid_min = -180. ;

lon:valid_max = 180. ;

lon:bounds = "lon_bounds" ;

lon:axis = "X" ;

int crs ;

crs:wkt = "GEOGCS[\"WGS 84\", \n  DATUM[\"World Geodetic System 1984\", \n    SPHEROID[\"WGS 84\", 6378137.0, 298.257223563, AUTHORITY[\"EPSG\",\"7030\"]], \n    AUTHORITY[\"EPSG\",\"6326\"]], \n  PRIMEM[\"Greenwich\", 0.0, AUTHORITY[\"EPSG\",\"8901\"]], \n  UNIT[\"degree\", 0.017453292519943295], \n  AXIS[\"Geodetic longitude\", EAST], \n  AXIS[\"Geodetic latitude\", NORTH], \n  AUTHORITY[\"EPSG\",\"4326\"]]" ;

crs:i2m = "0.002777777777778,0.0,0.0,-0.002777777777778,-180.0,90.0" ;

double lat_bounds(lat, bounds) ;

double lon_bounds(lon, bounds) ;

double time_bounds(time, bounds) ;

double time(time) ;

time:standard_name = "time" ;

time:long_name = "time" ;

time:axis = "T" ;

time:calendar = "standard" ;

time:units = "days since 1970-01-01 00:00:00" ;

time:bounds = "time_bounds" ;

// global attributes:

:id = "ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7cds" ;

:title = "Land Cover Map of ESA CCI brokered by CDS" ;

:summary = "This dataset characterizes the land cover of a particular year (see time_coverage). The land cover was derived from the analysis of satellite data time series of the full period." ;

:type = "ESACCI-LC-L4-LCCS-Map-300m-P1Y" ;

:project = "Climate Change Initiative - European Space Agency" ;

:references = "http://www.esa-landcover-cci.org/" ;

:institution = "UCLouvain" ;

:contact = "https://www.ecmwf.int/en/about/contact-us/get-support" ;

:comment = "" ;

:Conventions = "CF-1.6" ;

:standard_name_vocabulary = "NetCDF Climate and Forecast (CF) Standard Names version 21" ;

:keywords = "land cover classification,satellite,observation" ;

:keywords_vocabulary = "NASA Global Change Master Directory (GCMD) Science Keywords" ;

:license = "ESA CCI Data Policy: free and open access" ;

:naming_authority = "org.esa-cci" ;

:cdm_data_type = "grid" ;

:TileSize = "2025:2025" ;

:tracking_id = "e5dbaa6a-edbd-4b3e-b23f-b38055f4c4a5" ;

:product_version = "2.0.7cds" ;

:creation_date = "20181130T095525Z" ;

:creator_name = "UCLouvain" ;

:creator_url = "http://www.uclouvain.be/" ;

:creator_email = "landco...@uclouvain.be" ;

:source = "MERIS FR L1B version 5.05, MERIS RR L1B version 8.0, SPOT VGT P" ;

:history = "amorgos-4,0, lc-sdr-1.0, lc-sr-1.0, lc-classification-1.0,lc-user-tools-3.13,lc-user-tools-4.3" ;

:time_coverage_start = "20150101" ;

:time_coverage_end = "20151231" ;

:time_coverage_duration = "P1Y" ;

:time_coverage_resolution = "P1Y" ;

:geospatial_lat_min = "-90.0" ;

:geospatial_lat_max = "90.0" ;

:geospatial_lon_min = "-180" ;

:geospatial_lon_max = "180" ;

:spatial_resolution = "300m" ;

:geospatial_lat_units = "degrees_north" ;

:geospatial_lat_resolution = "0.002778" ;

:geospatial_lon_units = "degrees_east" ;

:geospatial_lon_resolution = "0.002778" ;

}

That were a lot of code blocks for my first question in this group. Any help is greatly appreciated!

Kind regards,

Bas

[Denis Sergeev]

Hi Bas,

Have you checked if those variables are stored as auxiliary coordinates of the main cube? I.e. what does print(cubelist[1]) give you?

Cheers,

Denis

[Klaus Zimmermann]

Hi Bas,

there are two issues with the file/iris and Denis already hinted at one.

Note that the file claims to follow CF conventions 1.6. I will
nevertheless refer to version 1.7; it is easy to find the corresponding
sections in the older document.

Really the whole file should end up in only one cube, coming from the
variable `lccs_class`. If you look at the ncdump output you will find
the attribute

```

lccs_class:ancillary_variables = "processed_flag current_pixel_state
observation_count change_count" ;

```
which tells us that those variables should be attached to the cube as
ancillary variables (see [1])
The bad news is that iris doesn't support this in the current release.
The good news is that one big PR has already been merged and remaining
work is ongoing so that iris 3 will support it.

Until that is released, perhaps the simplest way to fix this is to
remove the attribute so that all variables will be interpreted individually.


The other issue is with the `crs` variable. It's great that the data
provider gives us this crucial information about the underlying shape of
the earth! Unfortunately, it is not done in a CF compliant manner. The
details are described in [2,3,4]. If you are using something like
ESMValTool, you can fix this with a fix script. Basically you just have
to add an attribute to `lccs_class` (and the other variables if you
remove the ancillary connection) that says
```
lccs_class:grid_mapping = "crs" ;
```
and the crs should be modified to include
```
crs:grid_mapping_name = "latitude_longitude" ;
crs:crs_wkt = *
```
where * is the content of the current `crs:wkt` attribute.
In addition, the information from the well-known-text should be
replicated by appropriate attributes from [4] to give the full
information in a CF compliant manner.

Hope that helps!

Cheers
Klaus


[1]
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#ancillary-data

[2]
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#grid-mappings-and-projections

[3]
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#_latitude_longitude

[4]
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html,
Table F.1

On 01/11/2019 10:14, Bas Crezee wrote:
> Hi all, 
>
> I am using Iris to read landcover data recently released on the Climate
> Data Store. In case someone wants to reproduce the problem, I provide
> instructions for getting the data here. Below I copied the download
> script using the (Python) CDS api which I recommend, alternatively one
> can use the web form

> <https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-land-cover?tab=form> as
> well.
>
> |
> importcdsapi
> importos

>
>
> # NB: even though we request tar.gz, the downloaded file (as of 30 Oct
> 2019)
> # is netcdf and should be renamed to have extension .nc
>
>
>
>
> # Set parameters
> year =2015
> outdir ='/XXX/XXXX/landcover/'
>
>
> c =cdsapi.Client()
>
>

> filepath =os.path.join(outdir,'landcover_{0}.nc'.format(year))

>
>
> c.retrieve(
>     'satellite-land-cover',
>     {
>         'variable':'all',
>         'format':'tgz',
>         'version':'v2.0.7cds',
>         'year':str(year)
>     },
>     filepath)
>
> |
>
> I checked the file with the cfchecker

> <https://pypi.org/project/cfchecker/> and the format is largely OK. For

> Checkingvariable:time

> ------------------
>
>
> ERRORS detected:0
> WARNINGS given:1
> INFORMATION messages:2
>
> |
>

> *Now we come to my actual problem, it only reads in two variables
> properly, the others are simply neglected:*

> --
> You received this message because you are subscribed to the Google
> Groups "SciTools (iris, cartopy, cf_units, etc.) -
> https://github.com/scitools" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to scitools-iri...@googlegroups.com
> <mailto:scitools-iri...@googlegroups.com>.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/scitools-iris/d0e6ae0b-aec2-4d67-9808-51b051e9dfe0%40googlegroups.com
> <https://groups.google.com/d/msgid/scitools-iris/d0e6ae0b-aec2-4d67-9808-51b051e9dfe0%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Bas Crezee]

Hi Denis and Klaus,

Many thanks for your answer. Thanks in particular for your detailed explanation of the problem, Klaus! Your suggested workaround of removing that specific attribute works indeed. I also noted this messy layout of the `crs` variable, and I will report this back to the data provider.

Kind regards,

Bas

Op vrijdag 1 november 2019 11:57:04 UTC+1 schreef Klaus Zimmermann:

> an email to scitoo...@googlegroups.com
> <mailto:scitools-iris+unsub...@googlegroups.com>.

[Klaus Zimmermann]

Hi Bas,

great, thanks!

Cheers
Klaus

> <http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#_latitude_longitude>
>
>
> [4]
> http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html
> <http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html>,

> > CHECKING NetCDFFILE:landcover_2015.nc <http://landcover_2015.nc>

> <http://landcover_2015.nc>'

> > :creator_email = "landco...@uclouvain.be <javascript:>" ;

> > an email to scitoo...@googlegroups.com <javascript:>
> > <mailto:scitools-iri...@googlegroups.com <javascript:>>.

> > To view this discussion on the web, visit
> >
> https://groups.google.com/d/msgid/scitools-iris/d0e6ae0b-aec2-4d67-9808-51b051e9dfe0%40googlegroups.com
> <https://groups.google.com/d/msgid/scitools-iris/d0e6ae0b-aec2-4d67-9808-51b051e9dfe0%40googlegroups.com>
>
> >

> <https://groups.google.com/d/msgid/scitools-iris/d0e6ae0b-aec2-4d67-9808-51b051e9dfe0%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/scitools-iris/d0e6ae0b-aec2-4d67-9808-51b051e9dfe0%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
>
> --
> You received this message because you are subscribed to the Google
> Groups "SciTools (iris, cartopy, cf_units, etc.) -
> https://github.com/scitools" group.
> To unsubscribe from this group and stop receiving emails from it, send

> an email to scitools-iri...@googlegroups.com
> <mailto:scitools-iri...@googlegroups.com>.

> To view this discussion on the web, visit

> https://groups.google.com/d/msgid/scitools-iris/abb2ff41-350b-4281-8471-8ce98bfd3e4e%40googlegroups.com
> <https://groups.google.com/d/msgid/scitools-iris/abb2ff41-350b-4281-8471-8ce98bfd3e4e%40googlegroups.com?utm_medium=email&utm_source=footer>.