3

I successfully extracted some data from a netCDF file with a single variable E. (the original file was taking ages to load). I used the following code to extract the data:

ncks -d lon,24.,36. -d lat,21.,32. E_2015.nc extracted_E.nc

The resulting file I open using the popular python plugin "NetCDF Browser" by Mr Etienne Tourigny. Usually, the file opens and displays fantastically. However, I have come across some that have rather peculiar behaviour. (Please see below)

netCDF data in B/W Vector layer of where the data should be in purple

The purple layer is a vector layer for correct orientation reference, the netCDF data appears to be mirrored (with longitudes/latitudes switched and transversed).

gdalinfo gives me :

Driver: netCDF/Network Common Data Format
Files: extracted_E.nc
Size is 44, 48
Coordinate System is `'
Origin = (32.000000000000000,36.000000000000000)
Pixel Size = (-0.250000000000000,-0.250000000000000)
Metadata:
  E#long_name=Some Variable
  E#standard_name=Some Variable
  E#units=mm/day
  E#_FillValue=-999
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lat#_FillValue=-999
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  lon#_FillValue=-999
  NC_GLOBAL#history=Sun Nov 20 16:14:38 2016: ncks -d lon,24.,36. -d lat,21.,32. E_2015.nc extracted_E.nc
  NC_GLOBAL#NCO="4.5.4"
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={365,6}
  NETCDF_DIM_time_VALUES={16436,16437,16438,16439,16440,16441,16442,16443,16444,16445,16446,16447,16448,16449,16450,16451,16452,16453,16454,16455,16456,16457,16458,16459,16460,16461,16462,16463,16464,16465,16466,16467,16468,16469,16470,16471,16472,16473,16474,16475,16476,16477,16478,16479,16480,16481,16482,16483,16484,16485,16486,16487,16488,16489,16490,16491,16492,16493,16494,16495,16496,16497,16498,16499,16500,16501,16502,16503,16504,16505,16506,16507,16508,16509,16510,16511,16512,16513,16514,16515,16516,16517,16518,16519,16520,16521,16522,16523,16524,16525,16526,16527,16528,16529,16530,16531,16532,16533,16534,16535,16536,16537,16538,16539,16540,16541,16542,16543,16544,16545,16546,16547,16548,16549,16550,16551,16552,16553,16554,16555,16556,16557,16558,16559,16560,16561,16562,16563,16564,16565,16566,16567,16568,16569,16570,16571,16572,16573,16574,16575,16576,16577,16578,16579,16580,16581,16582,16583,16584,16585,16586,16587,16588,16589,16590,16591,16592,16593,16594,16595,16596,16597,16598,16599,16600,16601,16602,16603,16604,16605,16606,16607,16608,16609,16610,16611,16612,16613,16614,16615,16616,16617,16618,16619,16620,16621,16622,16623,16624,16625,16626,16627,16628,16629,16630,16631,16632,16633,16634,16635,16636,16637,16638,16639,16640,16641,16642,16643,16644,16645,16646,16647,16648,16649,16650,16651,16652,16653,16654,16655,16656,16657,16658,16659,16660,16661,16662,16663,16664,16665,16666,16667,16668,16669,16670,16671,16672,16673,16674,16675,16676,16677,16678,16679,16680,16681,16682,16683,16684,16685,16686,16687,16688,16689,16690,16691,16692,16693,16694,16695,16696,16697,16698,16699,16700,16701,16702,16703,16704,16705,16706,16707,16708,16709,16710,16711,16712,16713,16714,16715,16716,16717,16718,16719,16720,16721,16722,16723,16724,16725,16726,16727,16728,16729,16730,16731,16732,16733,16734,16735,16736,16737,16738,16739,16740,16741,16742,16743,16744,16745,16746,16747,16748,16749,16750,16751,16752,16753,16754,16755,16756,16757,16758,16759,16760,16761,16762,16763,16764,16765,16766,16767,16768,16769,16770,16771,16772,16773,16774,16775,16776,16777,16778,16779,16780,16781,16782,16783,16784,16785,16786,16787,16788,16789,16790,16791,16792,16793,16794,16795,16796,16797,16798,16799,16800}
  time#calender=standard
  time#long_name=time
  time#standard_name=time
  time#units=days since 1970-01-01 00:00:00 UTC
  time#_FillValue=-999
Corner Coordinates:
Upper Left  (  32.0000000,  36.0000000) 
Lower Left  (  32.0000000,  24.0000000) 
Upper Right (  21.0000000,  36.0000000) 
Lower Right (  21.0000000,  24.0000000)
(followed by band data)

Both the original netCDF and the extracted display correctly in netCDF browsers/viewers such as PanoplyJ. Other extracted netCDF files also display correctly in QGIS. Is there something special about this file? How can I fix this?

I have tried gdalwarp:

gdalwarp -t_srs '+proj=lonlat +datum=WGS84 +no_defs t_srs EPSG:4326' -of E_2015.nc x.nc

I get:

ERROR 1: Translating source or target SRS failed:

Furthurmore,

gdal_translate -a_ullr 24.0 32.0 36.0 21.0 -a_srs 'EPSG:4326' -of netCDF E_2015.nc x.nc

Gives:

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute

With the resulting image as shown below.

enter image description here

B/W file is the new data. This needs to be at right angles to what it is.

The resulting file when treated with the same code as above returns the original rotation (see first image above), but bound within the area of the purple vector reference. It is almost as if the Image is rotated around those two points the (ulx,uly) and (lrx,lry) overridden by gdal_translate (similar to flipping a card around it's diagonals). Also no amount of swapping x and y, or entering the opposing diagonal points in the gdal_translate command results in the correct orientation(on either the original data, or the result of gdal_translate).

I've also attempted to convert to Geotiff and reproject with no luck.

2

The netcdf does not have a CRS stored inside, so gdalwarp can not reproject to any CRS.

You better use

gdal_translate -a_srs 'EPSG:4326' -of netCDF E_2015.nc x.nc
  • Still the same orientation with output : Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute – mmfl Nov 20 '16 at 19:27
  • Can you share the file, or give a link to the original data? You might have to apply -a_ullr to set the file correctly. – AndreJ Nov 21 '16 at 6:31
  • Original files can be found here overriding with -a_ullr does change the orientation of the data by mirroring it across the x plane, but it still needs to be at right angles to the way it is now. Further translation with the resultant file doesn't adjust the data. is there a way to georeference the 4 corners rather than just ulx uly lrx lry? – mmfl Nov 21 '16 at 9:26
  • You can transform the file to tif and put it in the QGIS georeferencer. – AndreJ Nov 21 '16 at 20:14
  • Maybe related: gis.stackexchange.com/questions/185616/… – AndreJ Nov 21 '16 at 20:35
2

You can use the netCDF operators (NCO) to add spatial reference information that will be understood by QGIS. Essentially, we need to:

  1. Create a crs variable that has no associated dimensions
  2. Attach attributes to the crs variable that describe our coordinate system
  3. Use a grid_mapping attribute to associate the variables that contain spatial data to attributes defined under crs.

Here's an example in WGS84:

ncap -h -O -s 'crs=-9999' my_data.nc my_data.nc
ncatted -h -O \
    -a spatial_ref,crs,c,c,'GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.017453292519943295]]' \
    -a grid_mapping_name,crs,c,c,'latitude_longitude' \
    -a longitude_of_prime_meridian,crs,c,d,0 \
    -a semi_major_axis,crs,c,d,6378137 \
    -a inverse_flattening,crs,c,d,298.257223563 \
    -a grid_mapping,my_var_1,c,c,'crs' \
    -a grid_mapping,my_var_2,c,c,'crs' \
    my_data.nc
2

The problem is that the NetCDF file is not arranged in Climate and Forecast conventions and contains a grid of values without those being linked to a particular map projection and coordinates. However, looking at your gdalinfo output, it seems like the latitude and longitude values at each grid point are stored along with the grid of values.

A solution I worked out with similarly formatted NetCDF files:

1) Use a NetCDF utility to determine the lat/lon corner points. I've used ncdump for this though one of the NCO utilities should work as well, e.g.

ncdump -v latitude -f f (NetCDF file) | grep "latitude(1,maxy)"    # ULY
ncdump -v latitude -f f (NetCDF file) | grep "latitude(maxx,1)"    # LRY
ncdump -v latitude -f f (NetCDF file) | grep "latitude(1,1)"       # LLY
ncdump -v latitude -f f (NetCDF file) | grep "latitude(maxx,maxy)" # URY

repeat with longitude substituted above. maxx and maxy for your data is the maximum x and y values (44 and 48(?), respectively). You may not need all of these to georeference the data but might have a need later.

2) Once you find these coordinates, you'll need the upper left and lower right coordinates, i.e. longitude/latitude for 1,maxy and for maxx,1.

3) Knowing those, you can use gdal_translate to georeference the data in EPSG:4326 and save it in preferred format such as GeoTiff:

gdal_translate -a_srs EPSG:4326 -a_ullr ULX ULY LRX LRY (input) (output)

Performing gdalinfo on the resulting file will display the output file being georeferenced with the values used in gdal_translate.

4) You can use gdalwarp to then transform the resulting file into other projections.

NOTE: considering how the image seemed to have the x and y coordinates transposed when you used ncBrowse on the extracted file, you may have work with latitude and longitude values of the data. Just make sure since I can't grab the data you were working with.

Also, you may be able to avoid this if any of the NCO commands will also create the NetCDF file using CF conventions.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.