netcdf driver does not support rotated_pole projection

[hsumanto]

The netcdf driver does not recognized the rotated

Was: Fix the y-axis issue for projected CRS when translating from geotiff to netCDF (and vice versa from netCDF to geotiff)

This is a spin off ticket: #2893

Currently the driver (GDAL) expects top-down for projected crs but netCDF is bottom-up approach.

etienne has produced patch-geo2.txt in ticket #2129 which allow us to flip the y-axis, by applying the following around line 369 (in the part that deals with projected crs).

 /* netcdf standard is bottom-up, but leave it top first for now */
-       bBottomUp = FALSE;
+       bBottomUp = TRUE;

However, the resulting file will be upside-down in gdal and related software.

Thus, the total solution is required for projected crs so that the resulting file in geotiff -> netcdf conversion (and vice versa netcdf -> geotiff) will have the correct orientation for y-axis when viewed in both gdal and netcdf software.

[peifer]

Hi,

What is the current status of the y-axis issue? I am using gdal from trunk for reading netcdf files and all data come upside down. After reading the above, I changed line 716 in netcdfdataset.cpp as follows:

-    bBottomUp        = FALSE;
+    bBottomUp        = TRUE;

This helps me for the moment, but is hardly a sustainable "fix". Any real fixes planned for the foreseeable future?

PS: My netcdf files contain data from regional climate models. The data come in form of "rotated pole lat/lon grids", with the North Pole at 39.25N, 162W.

[etourigny]

We are in the midst of submitting an enhanced driver with full CF-1.5 Conventions export. Can you attach a subset of you file (less than 500kb) so we can make sure it works with your file?

[peifer]

My files are around 4M and come from an external data provider. To be honest, I wouldn't know how to reduce them to a suitable subset (without using GDAL). So I better leave them as they are. You can find 1 sample file here: ​http://transfer.eea.europa.eu/download/85bf0b19393459f66a3f658158a24961

[etourigny]

Dear peifer,

Thanks for reporting this issue, which is different from what this ticket is about (exporting from netcdf to geotiff), unless I don't understand how you are reading the data.

I have modified my dev version to assume bottom-up on import, I will notify you when you can test it.

Can you please tell me how you visualize the data? Are you converting to geotiff and if so, how? Which CRS are you assuming the data is in? I can't open the file with qgis as is (because it has 3 variables precip, lon and lat which GDAL sees as subdatasets).

I'd like to add that rotated pole grids are not yet supported by GDAL and PROJ.4 (which GDAL relies on), we are still working on that.

I downloaded the file, and it looks like it is not CF conformant. So I'm sorry to say that we cannot fully support this file. What is happening is that, since there is not valid coordinate system, GDAL just reads the file in the order in which the data was written. When we assume bottom-up by default, it isn't upside-down, but it's not well projected either (because of lack of support in GDAL for rotated grids, and the file is not CF compliant).

To show you why the file is not CF compliant, here is the output of 'ncdump -h TOT_PREC.lffd196101.dm.nc'

dimensions:
	x = 174 ;
	y = 190 ;
	time = UNLIMITED ; // (31 currently)
	rlat = 190 ;
	rlon = 174 ;
variables:
	double time(time) ;
		time:units = "day as %Y%m%d.%f" ;
		time:calendar = "proleptic_gregorian" ;
	float precip(time, y, x) ;
		precip:long_name = "total precipitation" ;
		precip:units = "mm/s" ;
	float lat(rlat, rlon) ;
		lat:standard_name = "latitude" ;
		lat:long_name = "latitude" ;
		lat:units = "degrees_north" ;
	float lon(rlat, rlon) ;
		lon:standard_name = "longitude" ;
		lon:long_name = "longitude" ;
		lon:units = "degrees_east" ;
	float rlat(rlat) ;
		rlat:axis = "Y" ;
		rlat:standard_name = "grid_latitude" ;
		rlat:long_name = "latitude in rotated-pole grid" ;
		rlat:units = "degrees" ;
	float rlon(rlon) ;
		rlon:axis = "X" ;
		rlon:standard_name = "grid_longitude" ;
		rlon:long_name = "longitude in rotated-pole grid" ;
		rlon:units = "degrees" ;

// global attributes:
		:CDI = "Climate Data Interface version 1.4.0.1" ;
		:Conventions = "CF-1.0" ;
		:history = "Fri Sep 10 09:25:36 2010: ncks -A latlon.nc Ptot_BCed_1961_1990_1961_1990_1961_01.nc\n",
			"Fri Sep 10 09:25:01 2010: ncks -v lat,lon /media/disk/DATA/ENSEMBLES_FLOODS/DMI-HIRHAM5//DMI-HIRHAM5_CTL_ERA40_FIX_25km_orog.nc ./latlon.nc" ;
		:Title = "total model precipitation" ;
		:CDO = "Climate Data Operators version 1.4.0.1 (http://www.mpimet.mpg.de/cdo)" ;
		:source = "HIRHAM" ;
		:experiment = "HC1" ;
		:institution = "DMI" ;
		:NCO = "20100910" ;
}

If you read the CF-1.0 specification at ​http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.0/cf-conventions.html#id2888877

a rotated pole CF file should look like:

dimensions:
  rlon = 128 ;
  rlat = 64 ;
  lev = 18 ;
variables:
  float T(lev,rlat,rlon) ;
    T:long_name = "temperature" ;
    T:units = "K" ;
    T:coordinates = "lon lat" ;
    T:grid_mapping = "rotated_pole" ;
  char rotated_pole
    rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
    rotated_pole:grid_north_pole_latitude = 32.5 ;
    rotated_pole:grid_north_pole_longitude = 170. ;
  float rlon(rlon) ;
    rlon:long_name = "longitude in rotated pole grid" ;
    rlon:units = "degrees" ;
    rlon:standard_name = "grid_longitude";
  float rlat(rlat) ;
    rlat:long_name = "latitude in rotated pole grid" ;
    rlat:units = "degrees" ;
    rlon:standard_name = "grid_latitude";
  float lev(lev) ;
    lev:long_name = "pressure level" ;
    lev:units = "hPa" ;
  float lon(rlat,rlon) ;
    lon:long_name = "longitude" ;
    lon:units = "degrees_east" ;
  float lat(rlat,rlon) ;
    lat:long_name = "latitude" ;
    lat:units = "degrees_north" ; 

The first problem is that your example file has the variable precip defined as float precip(time, y, x) ; and the axes x and y have no data associated with them, it should be defined as float precip(time, rlat,rlon) ;

It is also missing "grid_mapping" (rotated_pole) and "coordinates" information.

There is also an issue with the time axis which is not standard, but that doesn't have any effect in GDAL.

[peifer]

Thanks for looking so deeply into the matter and sorry for confusing issues. I thought this ticket would be appropriate as it mentions bBottomUp = TRUE; which was exactly what helped me as a brute-force quick fix for my "all data are upside down" issue.

Back to my netcdf: rotated pole grids ARE supported by PROJ.4 (through ob_tran transformation, in a somewhat counter-intuitive way, see my recent post at: ​http://osgeo-org.1803224.n2.nabble.com/forward-inverse-projection-logic-for-ob-tran-td6874351.html). So the CRS of the data is, in PROJ.4 notation:

+proj=ob_tran +o_proj=latlon +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180

However, GDAL/OGR's limited projection support via gdal/ogr/ogr_srs_proj4.cpp doesn't include ob_tran transformation, which is indeed a pity. I am hence forced to visualise and re-project my netcdf data in GRASS, which works fine. GRASS calls PROJ.4 for the re-projection and hands over the above parameters.

About the file not being CF conformant:
Thanks, I have taken note. I don't know much about netcdf standards and thought that the original data provider would know how to generate proper netcdf files. They don't, as it looks like :-(

[peifer]

What I forgot to mention: before importing the netcdf into GRASS, I generate a vrt file, in order to assign the correct bounding box and fix the mis-interpreted NODATA value (gdalinfo reports: NoData Value=9.96920996838686905e+36).

gdal_translate -of vrt netcdf:TOT_PREC.lffd196101.dm.nc:precip tmp.vrt \
               -a_ullr -22.27 21.01 16.01 -20.79 -a_nodata -9999

Then I bring the data into GRASS via r.in.gdal, which results into an upside-down view of Europe, which I then "fixed" with the bBottomUp = TRUE; trick.

[peifer]

Screenshots from GRASS for illustrating the issue

[etourigny]

Transferred this bug over to issues #4284 (import+export bottom-up by default) and #4285 (support for +proj=ob_tran within GDAL).

Closing this bug (marking as duplicate).