Avoid full slices

[ocefpaf]

pysgrid code perform full data slices in some places like:

https://github.com/sgrid/pysgrid/blob/master/pysgrid/sgrid.py#L497-L532

That can be memory expensive and/or slow when the data origin in remote.

Options:

1. Pass along the netCDF4-python variable instead of the slice. Maybe with little to no modification all the 2D processing will just work.
2. Make them a lazy array using either dask or biggus.

Option 2 adds a dependency and the smart thing to do would be to implement 1 having 2 as an optional behavior when the dependency is available. The advantages of 2 over 1 are: (a) to be able to load data that is bigger than memory; (b) avoid worrying about index logic when performing lazy computations.

[rsignell-usgs]

I wouldn't worry about the dependency and just go for dask (didn't the biggus guys say they were going to replace biggus with dask eventually?)

[ocefpaf]

There is a lot to do! This PR just avoid loading all the data in get_cell_center_lat_lon. However, pair_arrays will do that later 😒

We need either to make them lazy as soon as possible and adapt the code to use only lazy arrays, or to implement a custom __gettiem__.

[ocefpaf]

@ayan-usgs is the rational behind pairing the arrays? I know that pyugrid does the same, but I find that a violation of the data model. The netCDF file is not paired so why should we?

[rsignell-usgs]

If this came from pyugrid, perhaps @ChrisBarker-NOAA remembers the rationale...

[ocefpaf]

I asked him that at SciPy, but I was unconvinced 😉

I don't want to split hairs. And I do understand that it more compact to hold lon, lat into 1 variable.
However, we should strive for compatibility with the rest of the world.

Also, what is the point of doing this here to later on do this,

lon_data = grid_cell_centers[..., 0][sg_lon.center_slicing]
lat_data = grid_cell_centers[..., 1][sg_lat.center_slicing]

when plotting?

[rsignell-usgs]

I agree. I find that confusing. I have to then remember the order too.

[ocefpaf]

Maybe we just found another case for APIRUS 😉
@rsignell-usgs can you open an issue there for us please?

[ayan-usgs]

@ocefpaf,

There wasn't a strong motivating factor to do that other than to have some consistency with pyugrid at the time the code was written. I have no objection to changing that so that lat and lon are separate.

[ocefpaf]

Great @ayan-usgs! Let's coordinate to get this changed in both libraries then.

[ChrisBarker-NOAA]

BTW -- what are these pair_arrays you speak of?

[ocefpaf]

In pysgrid it is here:

pysgrid/pysgrid/utils.py

Lines 19 to 39 in 4f00871

	def pair_arrays(x_array, y_array):
	"""
	Given two arrays to equal dimensions,
	pair their values element-wise.
	For example given arrays [[1, 2], [3, 4]]
	and [[-1, -2], [-3, -4]], this function will
	return [[[1, -1], [2, -2]], [[3, -3], [4, -4]]].
	:param np.array x_array: a numpy array containing "x" coordinates
	:param np.array y_array: a numpy array containing "y" coordinates
	:return: array containing (x, y) arrays
	:rtype: np.array
	"""
	x_shape = x_array.shape
	paired_array_shape = x_shape + (2,)
	paired_array = np.empty(paired_array_shape, dtype=np.float64)
	paired_array[..., 0] = x_array[:]
	paired_array[..., 1] = y_array[:]
	return paired_array

[ChrisBarker-NOAA]

I see -- so this is essentially taking what may be two MxN arrays of, say, lat and lon, and making a MxNx2 array of lat-lon pairs.

OK -- I can see how that might start to feel a bit more ugly when you add extra dimensions...with 1-d it seems so natural.

Also -- a quick grep through the py_ugrid code does find a lot of "nodes[:,0]" etc....

But if you spilt lat and long up, you now could get them out of sync -- and do we make the rest of the API keep them apart as well, like the locate_nodes(self, points) method, etc.

Does anyone really like dealing with an individual point as a two distinct variables?

But the main driver here may be the lazy-loading thing -- if you want to lazy load, it's a lot easier to follow the data model of what you are loading from...

-CHB

[ocefpaf]

Does anyone really like dealing with an individual point as a two distinct variables?

I do. Also, xray, iris, virtually all plotting packages, and the CF-models is that way.

But the main driver here may be the lazy-loading thing.

Not really. That would be a bonus for this change. We could create a "lazy_paired" array.

it's a lot easier to follow the data model of what you are loading from

The motivation is not to make it easy, but consistent with the model we are loading from.

We have to remember that people will mix grid-like objects with data loaded with other packages (raw netCDF4-python, xray, iris, etc). IMO it is important to be consistent. Do we really want to see people plotting like?

ax0.plot(nodes[:,0], nodes[:,1], temp)
ax1.plot(ds.lon, ds.lat, ds.temp)