CDO

Hi Etienne!
Uwe is on a vacation. I guess, he will comment on this next week.

I have no personal experience with gdal, but it seems to be a widely used library.

• What are the use cases you have in mind?
• Apart from additional IO-Formats, what could be the benefit for CDO? You mentioned projections, but there seem to be a lot more.
• What is your concept of, for instance, fitting raster images into CDI's stream of gridded data?

Hi,

thanks for your answer, sorry for not responding sooner.

As a first comment, I would like to say that GDAL supports NetCDF and HDF so it is possible to convert raster data to NetCDF to use in cdo, and then back to another raster format (e.g. GeoTIFF) for further use in GIS software. However, I find it rather cumbersome.

Let me answer your questions:

1)
As GeoTIFF and other "raster" formats do not have a time axis, operators that work with temporal data would not be useful.

I was thinking of using cdo for its statistics, field stats and data manipulation operators (such as setc, setrtoc, add, addc, etc)
Recently I had to do simple comparisons in GeoTIFF files and ended up hacking the various existing GDAL python scripts to do it. I would have liked to use cdo instead.

2)
The raster formats supported by GDAL are listed here: http://www.gdal.org/formats_list.html

Thanks to the PROJ.4 library (http://trac.osgeo.org/proj/) GDAL supports:
- various geodetic datums (WGS84,NAD83,SAD69)
- various geographic (lat/lon) and projected (UTM, Albers, etc.) coordinate systems (referred to as CRS),
- transformations from one CRS to another
- various algorithms for CRS transformations and resampling

Therefore adding support for GDAL would enhance cdo's capacity for reprojecting and interpolation, as well as support for projected data. This would need some special consideration, as projected data is represented on a cartesian plane (x,y), versus geographic data which is in degrees.

GDAL also includes OGR (which deals with vector data (polygons), but I would not see any use for OGR in cdo.

3)
Raster data consists of binary data that can be represented as a matrix or array, so to fit them into a CDI stream should be "rather" straight-forward using the c api.
The grid definitions in CDI would have to be expanded to include the projected CRSs of PROJ.4 and GDAL.

The C and C++ apis provide raw RasterIO functions to read raster data into memory, similar to the netcdf library.
A simple example of raster access is at: http://www.gdal.org/gdal_tutorial.html

Raster files can have several Bands which I would see as various levels in the cdi stream. I do not have a great knowledge of CDI to comment further.

Regards,
Etienne

Hi Etienne,

We don't have experience with GDAL but as far as I know the support in CDO would be a very useful enhancement. I think its not possible to implement this support directly in CDI. Therefor I would propose to implement it as 2 CDO operator. E.g. import_gdal to read the GDAL files and export_gdal<,format> to write the CDO result with the GDAL library. <,format> is one of the GDAL supported raster formats, default could be geoTIFF. Here is an example for a processing step:

cdo export_gdal<,format>  [-op1 [-op2 [-opn]]] -import_gdal ifile ofile

Or simple convert a GDAL rasterfile to netCDF:

cdo -f nc import_gdal ifile ofile

Some examples of an import function are import_binary, import_cmsaf and import_amsr. An example for export is the undocumented function export_e5ml (import_e5ml).

You are right, one part of the work would be to expand CDI to support the projected CRSs of PROJ.4 and GDAL. CDI has limited support for projected CRSs (Sinusoidal, Lambert Azimuthal Equal Area and Lambert Conformal Conic). But the implementation is not very satisfacing. Therefor we plan to refactor this part of the code to have a common interface to describe projected grids. I think at this point we need the same interface for GDAL.

Is it possible to attach a few representiv example GDAL files? Particularly interesting is one with several Bands.

Best regards,
Uwe

PS: I'm on vacation also next week

Uwe,

I'm not sure how useful it would be to implement import and export to gdal, given that there is a gdal utility (gdal_translate) that does that, and gdal can support netcdf, hdf and grib.

My idea of supporting GDAL internally is to be able to work with various projected CRSs (such as UTM), without having to convert to and from geographic coordinates as well as to/from netcdf files.

However, these import/export operators would make it easier for who does not want to learn GDAL, and provide support for the types that GDAL does not support (??) as well as some that cdo might not (e.g. HDF files from NASA if I'm not mistaken).

As long as you guys are planning to refactor the treatment of projected grids, perhaps it would be a good opportunity to use PROJ.4 and/or GDAL internally? GDAL has useful functions for querying the CRS info, such as IsProjected()

I've attached a zip with a few geotiff files, and one netcdf file (converted with gdal_translate).
- TM.burnpix.PNSCa.2009.tif TM.burnpix.PNSCa.2010.tif TM.burnpix.PNSCa.2009-10.tif
the last has 2 bands (created from the 2 other files), all in UTM coordinates
- TM.burnpix.PNSCa.2009-10.wgs84.tif in lat/lon coordinates
- TM.burnpix.PNSCa.2009-10.nc and TM.burnpix.PNSCa.2009-10.wgs84.nc
converted to netcdf, has an issues with lat/lon specifications not adhering to COARDS though (trying to get the fix in the official release)

PS: don't trouble yourself with this too much for now, you're on vacation!

Hi Etienne,

Thanks for the geotiff files! I found some time to play a little with GDAL. The biggest problem for the implementation of the internal support (CDI) would be to automaticaly detect the file format and choose the corresponding I/O library. E.g. a netCDF file could be read either directly with the netCDF library over CDI or with GDAL?

I will start with the refactoring of the support for projected CRSs in June. The current plan is to extend the grid type curvilinear with the following functions:

int gridIsProjected(gridID);
void gridDefProjName(int gridID, char *projectionName);
void gridInqProjName(int gridID, char *projectionName);
void gridInqProjParNum(int gridID);
void gridDefProjPar(int gridID, char *parName, double param);
void gridInqProjPar(int gridID, int index, char *parName, double *param);

So it would be possible to define the name and all parameters for a projection. The transformation to geographical coordinates is not part of the I/O library. We will do this transformation with PROJ.4 directly in CDO.

Best regards,
Uwe

Uwe -

sorry for the late reply. My advice with netCDF would be to use internal libnetcdf. The GDAL netCDF driver uses libnetcdf and does some stuff with coordinates which would lead to inefficiencies and unexpected results (when used inside libcdi/cdo).

If I understand correctly, you want to define all projected CRSs as type GRID_CURVILINEAR? The following CDO grids are supported by PROJ.4/GDAL: GRID_LCC, GRID_LCC2, GRID_LAEA, GRID_SINUSOIDAL .

GDAL allows to specify a CRS using complete WKT, PROJ.4 or EPSG:n format, you might want to look at that.
http://www.gdal.org/gdalinfo.html
http://www.gdal.org/ogr/osr_tutorial.html
http://www.gdal.org/gdal_tutorial.html

It makes sense to start with PROJ.4 as it is simpler to install. I would recommended adding libgeotiff to the I/O library. As geotiffs hold CRS and geotransform information, it would be sufficient in my view. Therefore, passing to and from other raster formats supported by GDAL would be a breeze.

Perhaps a final step would be the integration with GDAL as to incorporate functionality not in PROJ.4 (e.g. more options for CRS definition, file formats, remapping methods).

regards, Etienne

The function OGRSpatialReference::SetFromUserInput(const char * pszDefinition) from GDAL/OGR handles various user input

http://www.gdal.org/ogr/classOGRSpatialReference.html#ec3c6a49533fe457ddc763d699ff8796

Uwe,

Any chance you had time to look into your proposed support for projected CRS?

Here is how I think it should work

- support PROJ.4 (which also supports EPSG codes) via projectionName variable as you suggested
- add GeoTIFF reading and writing support to CDI with the libgeotiff library
- any other formats should be converted to/from GeoTIFF with GDAL
- optional: WKT can be translated to/from PROJ.4 strings outside of cdo by the user, or with a helper script (using OGR from GDAL)

If there is any way I can help, please let me know!

Etienne

Hi Uwe and Etienne,

thank you very much for this chat!
It seems to be long time ago you have worked on it - but recently I am trying to solve something where I need something related to this. And if there would be a "CDO" solution I would be extremely happy.

So my question is related to this CRS issue, as I need to get an EPSG code description into a netcdf file. In other words, "project" a geo-lat/lon grid in a Lambert Equal Area information already in the file, which is needed for a GEOSERVER.

The nc-header needs someething like this:

netcdf WindSpeed_Current_t {                                                                                
dimensions:                                                                                                 
        lon = 344 ;                                                                                         
        lat = 168 ;                                                                                         
variables:                                                                                                  
        float Band1(lat, lon) ;                                                                             
                Band1:long_name = "GDAL Band Number 1" ;                                                    
                Band1:_FillValue = -9999.f ;                                                                
                Band1:grid_mapping = "crs" ;                                                                
        char crs ;                                                                                          
                crs:grid_mapping_name = "latitude_longitude" ;                                              
                crs:longitude_of_prime_meridian = 0. ;                                                      
                crs:semi_major_axis = 6378137. ;                                                            
                crs:inverse_flattening = 298.257223563 ;                                                    
                crs:spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.2572
23563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",
\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" ;
                crs:GeoTransform = "-38.24372055771543 0.2847632482918333 0 72.0355784629853 0 -0.2839812015
549089 " ;

Could/would this "cdo export_gdal<,format> [-op1 [-op2 [-opn]]] -import_gdal ifile ofile" command do so?
What are the options? And would it work for importing "normal" netcdf?

Sincerely Yours
Swantje