Writing netcdf from python changes pixel/gridline node registration

Q&A PyGMT Q&A
1

How can we write pixel-node-registration netcdf files in Python? Using a python netCDF library, I’m creating files that claim to be pixel-node registration, but they’re actually line node registration. It’s causing the grid to become warped by 1/2 pixel on each side every time I read/write. I can’t find any mention of node registrations in third party NetCDF libraries. Any help would be greatly appreciated!

In the example below, I read and write a simple .grd file.

  • In GMT, original_netcdf.grd file is pixel node registration. gmtinfo -R correctly gives range limits at the outside edges of the leftmost and rightmost pixels.
  • The python library reads it as pixel node, with one array element per pixel. The edges of the python array are offset by half a pixel inwards from gmtinfo -R, since they correspond to the centers of each pixel. The dimensions of the array in Python are correct.
  • Finally, the python library writes the array wrong. When it writes to netcdf, it assumes that the first and last elements of the array correspond to the edges of gmtinfo -R, as if it were gridline node registration, when in reality they correspond to the pixel centers. This squeezes the entire array in real space by about half a pixel, but keeps the array dimensions the same (which is why I didn’t notice this bug for years!). Every pixel gets ever so slightly smaller and the increment gets ever so slightly smaller, as can be seen in the gmtinfo calls below.

Note: calling grdsample -T would of course cause the dimensions to change by 1, which I don’t want.

Setup

kmaterna-MBP-3:$ gmt grdinfo original_netcdf.grd
original_netcdf.grd: Pixel node registration used [Geographic grid]
original_netcdf.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
original_netcdf.grd: **x_min: 242.55 x_max: 246 x_inc: 0.00166666666667** name: longitude [degrees_east] n_columns: 2070
original_netcdf.grd: y_min: 31.9166666667 y_max: 35.8888888889 y_inc: 0.00138888888889 name: latitude [degrees_north] n_rows: 2860
original_netcdf.grd: z_min: 0.00355224008672 z_max: 0.984841525555 name: z
original_netcdf.grd: scale_factor: 1 add_offset: 0
original_netcdf.grd: format: netCDF-4 chunk_size: 130,130 shuffle: on deflation_level: 3

Action

Then in Python, I do:

[x, y, z] = read_netcdf4_python('original_netcdf.grd');
write_netcdf4_python(x, y, z, 'python_written_file.nc');

Result

kmaterna-MBP-3:$ gmt grdinfo python_written_file.nc 
python_written_file.nc: Pixel node registration used [Cartesian grid]
python_written_file.nc: Grid file format: nd = GMT netCDF format (64-bit float), COARDS, CF-1.5
python_written_file.nc: **x_min: 242.550833333 x_max: 245.999166667 x_inc: 0.00166586151369** name: x [range] n_columns: 2070
python_written_file.nc: y_min: 31.9166666665 y_max: 35.8888888887 y_inc: 0.00138888888888 name: y [azimuth] n_rows: 2860
python_written_file.nc: z_min: 0 z_max: 0 name: z [mm]
python_written_file.nc: scale_factor: 1 add_offset: 0
python_written_file.nc: format: classic

Python Specs:

I’m using the netCDF4 python library with Python 3.7. The problem also exists with the scipy.io.netcdf library. I can provide the text of the read/write functions in python if that would help.

2

update: this might be library/dependency problems beyond my pay grade. The registration/range/increment of the output file depends on whether Python writes as float32 vs float64, AND it depends on which computer I’m using. I’ve created at least 4 different outcomes by testing these variables. So, I’m likely to use the combination that works for my needs and not poke the bear any further.

The most annoying problem is that from Python, writing the file as float32 consistently results in gridline node registration while writing the file as float64 results in pixel node registration. Weird.

3

The notion of pixel vs grid node registration has nothing to do with netCDF. For example in GeoTTif world they are known as pixel-is-area and pixel-is-point.

Also you are using functions read|write_netcdf4_python that we know nothing about.

4

Thank you for letting me know. I asked here because I couldn’t find any equivalent of pixel vs grid for general netcdf files… now I understand why. The behavior I’m seeing must be an unexpected quirk of the python/netcdf interface.

5

People are always welcome to ask things here … as long as we are supposed to be able to help. But using functions that we (well, I) know nothing about doesn’t fit that bill.

6

How about changing the NetCDF header using grdedit -T to toggle between gridline and pixel registration?

Or, you could take a look at salem. it’s an xarray accessor library that has a salem.Grid.corner_grid and salem.Grid.center_grid and you can convert between the two. See also this PyGMT issue on Deciding on default grid registration and grid type for PyGMT · Issue #487 · GenericMappingTools/pygmt · GitHub that talks a bit about the Pixel/Gridline problem. Most Python libraries I’m familiar with that handle NetCDF assume pixel registration, as they’re usually written by optical remote sensing/earth observation scientists :slightly_smiling_face:.