transpose large netCDF with GDAL, NCL, NCO or CDO

Question

I have a large netCDF file (444GB) that I would like to transpose. The file is flipped with latitude on the x-axis and lon on the y-axis. See pictures of one band below.

Not correct

Correct: I transposed this using raster package in R but my 444GB file is too big and the process has been running for 30 hours now non-stop. Thus, it is too slow to accomplish this task. How can I achieve the same task using GDAL, NCL, NCO, or CDO?

An example file can be gotten from here link to files

Answer 1

You mention using R - would a Python solution be suitable? I'm not sure how this will compare in terms of processing time.

The netCDF4 Python library allows you to store your variable in a numpy array. You can then use numpy to manipulate the array, in your case you can use the numpy.transpose() method.

The following shows how to access a 2D array from a netCDF file and save it as a GeoTIFF file:

import gdal,gdalconst,osr
import numpy
from netCDF4 import Dataset

# Open the NetCDF file for reading.
fh = Dataset(r"INPUT_PATH.nc","r")

# Read your parameter into a numpy array.
# Depending on your dataset you may have to specify additional index values
# e.g. Time.
values = numpy.array(fh.variables["YOUR_PARAMETER"])

# Now you can manipulate your array.
"""
Transpose method goes here
"""

# Get array dimensions so we can create an equal size tiff file.
rows,cols = values.shape

# Create output tiff file.
# driver.Create() parameters are: output path, number of columns, number of rows,
# number of bands, data type
driver = gdal.GetDriverByName("GTiff")
out_data = r"OUTPUT_PATH.tif"
out_tif = driver.Create(out_data, cols, rows, 1, gdalconst.GDT_Float32)

# Create Spatial Reference object and set GeoTIFF projection.
# This information may be found in either the data documentation or the netCDF file.
prj = osr.SpatialReference()
prj.ImportFromEPSG(4326) # WGS84
out_tif.SetProjection(prj.ExportToWkt())

# Set GeoTransformation.
# This information may be found in either the data documentation, the netCDF file, or
# can be derived. For example, if you know the longitude range and number of columns
# you can calculate the x step as float(lon_range)/float(num_cols).
geotrans = [top_left_x,x_step,rotation,top_left_y,rotation,y_step]
out_tif.SetGeoTransform(geotrans)

# Finally we can write the array to the raster band.
out_band = out_tif.GetRasterBand(1)
out_band.WriteArray(values)

# Clear the memory and close the output file.
out_tif.FlushCache()
out_tif = None

Your main analysis will go in the "Transpose method goes here" section. When transposing the array remember that you will need to account for these changes in the output GeoTIFF creation and geotransformation. For example the number of rows and columns will be reversed, as will the starting coordinates and x and y steps.

Answer 2

In NCO, the ncpdq operator can transpose (and flip) dimensions as needed:

ncpdq --rdr=dim2,dim1 incorrect_file.nc correct.nc

The relevant section of the documentation at http://nco.sourceforge.net/nco.html#ncpdq-netCDF-Permute-Dimensions-Quickly is:

ncpdq re-shapes variables in input-file by re-ordering and/or reversing dimensions specified in the dimension list. The dimension list is a whitespace-free, comma separated list of dimension names, optionally prefixed by negative signs, that follows the ‘-a’ (or long options ‘--arrange’, ‘--permute’, ‘--re-order’, or ‘--rdr’) switch. To re-order variables by a subset of their dimensions, specify these dimensions in a comma-separated list following ‘-a’, e.g., ‘-a lon,lat’. To reverse a dimension, prefix its name with a negative sign in the dimension list, e.g., ‘-a -lat’. Re-ordering and reversal may be performed simultaneously, e.g., ‘-a lon,-lat,time,-lev’. ... To create the mathematical transpose of a variable, place all its dimensions in the dimension re-order list in reversed order. This example creates the transpose of three_dmn_var:

% ncpdq -a lon,lev,lat -v three_dmn_var in.nc out.nc
% ncks --trd -C -v three_dmn_var in.nc
lat[0]=-90 lev[0]=100 lon[0]=0 three_dmn_var[0]=0 
lat[0]=-90 lev[0]=100 lon[1]=90 three_dmn_var[1]=1 
lat[0]=-90 lev[0]=100 lon[2]=180 three_dmn_var[2]=2 
...
lat[1]=90 lev[2]=1000 lon[1]=90 three_dmn_var[21]=21 
lat[1]=90 lev[2]=1000 lon[2]=180 three_dmn_var[22]=22 
lat[1]=90 lev[2]=1000 lon[3]=270 three_dmn_var[23]=23 
% ncks --trd -C -v three_dmn_var out.nc
lon[0]=0 lev[0]=100 lat[0]=-90 three_dmn_var[0]=0
lon[0]=0 lev[0]=100 lat[1]=90 three_dmn_var[1]=12
lon[0]=0 lev[1]=500 lat[0]=-90 three_dmn_var[2]=4
...
lon[3]=270 lev[1]=500 lat[1]=90 three_dmn_var[21]=19
lon[3]=270 lev[2]=1000 lat[0]=-90 three_dmn_var[22]=11
lon[3]=270 lev[2]=1000 lat[1]=90 three_dmn_var[23]=23