1D "GRID" issues

[bjarnebuchmann]

Dear Wizard(s) (@tsutterley),

I am trying to use pyTMD to create boundary conditions for a hydrodynamic model. Therefore, I need tidal elevations and depth-integrated currents in a number of points - and for a number of time steps. Typically a few hundred points times a few hundred time values.

I would prefer to use the "scripts"-approach, ie compute_tidal_elevations.py and compute_tidal_elevations.py.

Presently, it seems, however, that it is assumed that the target should be (ny,nx,nt)-shaped, ie that input data are always given as 1D coordinates. However, I am really interested in x,y (lon.lat) being defined in specific points (along a track/boundary), and thus they should not be "crossed" to a "grid". Does this make sense?
In essence, I give lon(bdy), lat(bdy), time(time), and I would like output to be shaped (bdy,time), ie 2D.

It would not be a problem for me to create a fake/degenerate dimension of length 1, but I cannot position my points to be on a proj4 grid.

Any ideas on how to approach this?
Would I have to write my own wrapper-script?

Thank you,
/Bjarne

[bjarnebuchmann]

This is not an answer, but rather that an answer is not really that the 1D common dimension does not seem possible with the present implentation of spatial.py and compute_tidal_elevations.py

The closest to my use case seems to be compute_tidal_{elevations,currents}.py with --format=netCDF4 --type=grid.
However, compute_tidal_elevations.py does not look at the dimensions of input x, y, but defines them as separate coordinates variables, even if I specify input with a common single dimension for x and y.
For instance, if I specify

dimensions:
        time = 7 ;
        ibdy = 10;
variables:
        float lon(ibdy) ;
                lon:_FillValue = -9999.f ;
                lon:long_name = "longitude at T-points (cell centres)" ;
                lon:units = "degrees-east" ;
        float lat(ibdy) ;
                lat:_FillValue = -9999.f ;
                lat:long_name = "latitude at T-points (cell centres)" ;
                lat:units = "degrees-north" ;
.. ..

Then the output data becomes

dimensions:
	y = 10 ;
	x = 10 ;
	time = 7 ;
variables:
	float time(time) ;
		time:long_name = "Time" ;
		time:units = "days since 1992-01-01T00:00:00" ;
		time:calendar = "standard" ;
	double lon(y, x) ;
		lon:long_name = "Longitude" ;
		lon:units = "Degrees_East" ;
	double lat(y, x) ;
		lat:long_name = "Latitude" ;
		lat:units = "Degrees_North" ;

The issue here lies in that the data are read without checking dimensionality (spatial.py, line 345):

 dinput[key] = fileID.variables[nc][:]

So, data are read as if they are 1D, and without checking if they share dimension.
It is thus assumed that they have separate dimensions, and subsequently, they are passed to meshgrid(x,y) (compute_tidal_elevations.py , line 270).

To solve my use case, it would be necessary to look at the actual dimensions of x and y, and only do a meshgrid if the dimensions do not match. Thus, I presently believe that I cannot use the pyTMD scripts to solve my use case. I might be able to write a different/specific wrapper script. But I would have to also re-code how netCDF files are read, because the dimension-info is not retained in from_netCDF4().

By the way @tsutterley: Is there a particular reason for having the attributes_list and not simply grab and copy all attributes from the input file variables? fileID.variables[nc].__dict__ should give you all attributes as a dict, ready to go. If you want, you can the n use if attr in dict, which seems simpler than the try:/except: case(?)

Best,
/Bjarne

[tsutterley]

I'm adding functionality for the "time series" type in #162. This has also been an outstanding, undocumented issue that's been on the list. Hopefully this works for your use case.
No particular reason for not getting all of the attributes. In other programs I do iterate over ncattrs(). Might change it in the future to be similar.
-ts

[bjarnebuchmann]

Great.
I'll look forward to testing it when ready (So far, still just grid or drift options for --type.)

[bjarnebuchmann]

I've recently started to use the following "trick" for passing var:attrs from src to tgt nc files. It is important to treat _FillValue special, thus the pop. Note that the scheme works even if the src var does not have a _FillValue attribute:

# Get attrs from src var:
vattrs=ncsrc.variables[vnam].__dict__
vtyp=ncsrc.variables[vnam].dtype
# Get fillvalue if it is here, otherwise get None
fillval=vattrs.pop('_FillValue',None)
# Note that fillval may be None in the following:
tgtvar=nctgt.createVariable(vnam,datatype=vtyp,dimensions=WhateverDimList,fill_value=fillval)
# Transfer all the rest of the attributes (excepting fillval, which was popped above):
tgtvar.setncatts(vattrs)
# nctgt.variables[vnam].setncatts(vattrs) # Same, if tgtvar is not stored