Converting NOAA's World Ocean Database NetCDF files to tabular data using Python

Question

I'm trying to convert variables from NetCDF files to tabular data, but can't understand how to group observations by some index or key variable.

For example, when loading data from WOD queries using the xarray library

import xarray as xr

ds = xr.open_dataset(r'ocldb1573479266.24760_OSD.nc')
ds

I get something like

<xarray.Dataset>
Dimensions:                         (Alkalinity_obs: 32749, Ammonia_obs: 15252, CFC11_obs: 16256, CFC12_obs: 16317, Chlorophyll_obs: 58045, DeltaC14_obs: 1456, DeltaHe3_obs: 1929, Helium_obs: 2143, Neon_obs: 377, Nitrate_obs: 108954, Oxy18_obs: 409, Oxygen_obs: 300929, Phosphate_obs: 183091, Pressure_obs: 581435, Salinity_obs: 581435, Silicate_obs: 167977, Temperature_obs: 581435, Tritium_obs: 2015, biosets: 1149, casts: 51401, numberofpis: 6079, pCO2_obs: 8513, pH_obs: 42365, tCO2_obs: 24672, z_obs: 581435)
Coordinates:
    lat                             (casts) float32 ...
    lon                             (casts) float32 ...
    time                            (casts) datetime64[ns] 1890-06-27T02:00:00.000214528 ... 1988-03-05
    z                               (z_obs) float32 ...
Dimensions without coordinates: Alkalinity_obs, Ammonia_obs, CFC11_obs, CFC12_obs, Chlorophyll_obs, DeltaC14_obs, DeltaHe3_obs, Helium_obs, Neon_obs, Nitrate_obs, Oxy18_obs, Oxygen_obs, Phosphate_obs, Pressure_obs, Salinity_obs, Silicate_obs, Temperature_obs, Tritium_obs, biosets, casts, numberofpis, pCO2_obs, pH_obs, tCO2_obs, z_obs
Data variables:
    country                         (casts) |S40 ...
    WOD_cruise_identifier           (casts) |S40 ...

(. . .)

    Temperature                     (Temperature_obs) float32 ...
    Temperature_sigfigs             (Temperature_obs) int8 ...
    Temperature_row_size            (casts) float64 7.0 6.0 7.0 ... 58.0 45.0
    Temperature_WODflag             (Temperature_obs) int8 ...
    Temperature_origflag            (Temperature_obs) float32 ...
    Temperature_WODprofileflag      (casts) int8 ...
    Temperature_Access_no           (casts) float64 ...
    Temperature_Scale               (casts) |S170 ...
    Temperature_Instrument          (casts) |S170 b'' b'' b'' ... b'' b''

(. . .)

When inspecting shapes

print('CASTS          : ', ds.casts.to_dataframe().shape )
print('Temperature    : ', ds.Temperature.to_dataframe().shape )
print('Temperature_obs: ', ds.Temperature_obs.to_dataframe().shape )
print('time           : ', ds.time.to_dataframe().shape )

I get the following

CASTS          :  (51401, 4)
Temperature    :  (581435, 1)
Temperature_obs:  (581435, 1)
time           :  (51401, 4)

So, as expected, it seems for each cast there is a set of Temperature observations. But I'm not familiar with xarray (or NetCDF4-python) methods for mask, group, merge, etc. So my question is: how can I create something like this?

|     | cast | lat | lon | temperature | ... |
|:---:|:----:|:---:|:---:|:-----------:|:---:|
|  0  |   0  |     |     |             |     |
|  1  |   0  |     |     |             |     |
|  2  |   0  |     |     |             |     |
| ... |      |     |     |             |     |
|  21 |   0  |     |     |             |     |
|  22 |   1  |     |     |             |     |
|  23 |   1  |     |     |             |     |

I already tried to use to_dataframe() method, but it did not work.

Answer 1

Since there is none answer since question was posted, it might be helpful to tell how I dealt with this issue.

First of all, I'm still not familiar with .netcdf format and therefore it took me some time to suspect that the problem wasn't exclusively this lack of familiarity nor some intrinsic complexity of the format, but the apparent absence of an explicit key/id variable to merge stuff. It is quite clear that this variable it should be casts, but I couldn't find a way to use it in order to create flat tables. Additionally, I made tests with netcdfs from other sources and it was trivial to manipulate information because there were "id" variables in them. Then I suspect there is an issue related to the way data come from WOD requests.

That said, what solved my problem was wodpy:

# i am sorry for the aligned signs in advance, it is a bad practice, but an old habit.

from wodpy import wod

FILES   = [
'ocldb1573755495.4900.OSD'
'ocldb1573831711.15486.OSD'
'ocldb1573836294.18286.OSD'
]


for fname in FILES:
  # load files from download from wod 
  fid = open("{}".format(fname), newline=None)
  p   = wod.WodProfile(fid)
  p.return_file_position_to_start_of_profile(fid)

  # create dataframe
  _df        = p.df().copy()
  _df['uid'] = p.uid() # thats the id
  """
  do stuff here
  """
  df = _df.copy()

  while p.is_last_profile_in_file(fid) == False:
    try:
      del _df
      p          = wod.WodProfile(fid)
      _df        = p.df().copy()
      _df['uid'] = p.uid()
      """
      do stuff here
      """
      df         = df.append(_df)

    except:
      pass

Probably there are clever ways to do this, mainly because hundreds of errors arise during the loop above, then data is fatally lost. Nevertheless it was good enough for my needs at the occasion.

Lastly, even with this approach something went wrong at specific point of loop. My requests were measurements from 1950 to 2010, file chunks came sorted by date, then it was possible to find out that something in 2005 interrupted the loop in such a way that no error was raised. Different files from different requests of same time range froze the loop at same point. So it's probably a thing to take a look in case of apparent infinite loop.