xgrid.max() and xgrid.min() are returning the same value in readers/interpolation/interpolators.py

[arian-dialectaquiz]

Hello. With your guidelines to my previous issue on pixelsize I was able to fix the problem by creating code that imports the modelgrid from the Eulerian model and excludes NaNs. And the reader created for the ECOM (uploaded .txt file) model worked together with the others.

I recently had to reinstall Opendrift and I ran into the following problem:

In the Linear2DInterpolator () class in readers / interpolation / interpolators.py there is a division by zero in the definitions of self.xi and self.yi, since maximum and minimum xgrid and ygrid are returning as having the same values, what makes all subsequent values become NaN in the interpolation (first figure).

I did some tests, returning to the previous mode of the reader, ie, without the codes that exclude NaNs from the model grid, and returning the division by pixelsize in basereader / variables.py and ended up facing the same problem as xgrid.min () == xgrid.max () and ygrid.min == ygrid.max ().

Could you give me some guidance on what is possibly causing this problem that did not occur a week ago? Could it be possible that I installed something wrong when reinstalling Opendrift?

Thank you in advance. reader_ECOM_S2Z.txt

[knutfrode]

Hi,

Since you do not have/know any projection here, I presume the lon and lat arrays are 2D? And as mentioned, these have to be valid (and continuous) for all points, they cannot be NaN e.g. over land.

In this case, you should not define the xmin, xmax, delta_x (same for y), as these are now determined automatically by the parent class (previously it was done in e.g. ROMS reader). You are also defining xgridand ygrid, but it does not seem like these are used for anything?

Thus I believe you may remove the definitions of these mentioned attributes of the self object. Does it work if you remove this?

[arian-dialectaquiz]

Hello. The xgrid and ygrid definitions were my attempt to create an array going from xmin to xmax, while delta_x, so that in interpolators.py the functions .min () and .max () would not return the same result.

However the problem remains the same after excluding the definitions of xmax, xmin, delta_x and xgrid (the same for y).

What are xgrid and ygrid, I did not find their definitions in the interpolations readers, not in the basereaders. And, no other reader (like ROMS native) actually declares them as you sad.

[knutfrode]

I would then make sure to remove everything regarding xmax, xmin, delta_y and x_grid (same for y) from your reader, to avoid any more confusion from these.

So the only thing that matters here are the arrays self.lon and self.lat Are these continuous and without holes/NaN? Do you get reasonable output about geographical coverage from print(reader)(before starting the simulation) ? Could you post this output?

[arian-dialectaquiz]

I removed everything regarding the xamx, xmin, delta_x and xgrid (same for y), mantaining just the self.lat and self.lon and indx = np.floor(x).astype(int) + clipped indy = np.floor(y).astype(int) + clipped.

The geographial coverage output from print(reader) is:

Where lat and lon are 2D arrays continuous and wihtout NaN.

There are NaN in depth values, but they are removed later in de reader. This might be causing the problem?

However the problem keeps showing.

[knutfrode]

The output from >>> print(reader) should look something like this:

===========================
Reader: roms native
Projection: 
  None
Coverage: [pixels]
  xmin: 0.000000   xmax: 150.000000   step: 1   numx: 151
  ymin: 0.000000   ymax: 80.000000   step: 1   numy: 81
  Corners (lon, lat):
    (  5.41,  67.56)  ( 16.39,  71.68)
    ( 11.10,  65.64)  ( 22.02,  69.34)
Vertical levels [sigma]: 
  [-0.98571429 -0.95714286 -0.92857143 -0.9        -0.87142857 -0.84285714
 -0.81428571 -0.78571429 -0.75714286 -0.72857143 -0.7        -0.67142857
 -0.64285714 -0.61428571 -0.58571429 -0.55714286 -0.52857143 -0.5
 -0.47142857 -0.44285714 -0.41428571 -0.38571429 -0.35714286 -0.32857143
 -0.3        -0.27142857 -0.24285714 -0.21428571 -0.18571429 -0.15714286
 -0.12857143 -0.1        -0.07142857 -0.04285714 -0.01428571]
Available time range:
  start: 2016-02-02 12:00:00   end: 2016-02-04 12:00:00   step: 1 day, 0:00:00
    3 times (0 missing)
Variables:
  sea_ice_area_fraction
  sea_floor_depth_below_sea_level
  sea_ice_thickness
  land_binary_mask
  land_binary_mask
  sea_water_salinity
  sea_water_temperature
  x_sea_water_velocity
  sea_ice_x_velocity
  y_sea_water_velocity
  sea_ice_y_velocity
  sea_surface_height

So probably you are printing something else(?)

The only problem I can identify in your log is that a data-block of size 1x1x1 is fetched, whereas this should be at least 5x5 pixels large. So there must be something wrong also in the get_variables() method in your reader.

[arian-dialectaquiz]

Yes, I sorry. I printed my own report.

Now the print(reader):

I can change the data-block by changing the bufer size, right? I'll look to the get_variables() method, comparing it to the ROMS native one.

[knutfrode]

Yes, this looks quite ok.

The times and time step are a bit strange (2 hours 28 minutes), but this should not be a problem in itself.

There are also some variables listed, which are more internal variables, and should probably be removed, such as e.g. "X_coordinate_in_Cartesian system, corner_longitude and angle_of_rotation_from_east_to_x. These are more internal parameters, whereas the reader should rather return variables that are needed for drift calculations (current, wind, waves etc).

Anyway this indicates that your __ init__ method is probably ok, and that the problems might be in the get_variablesmethod. So yes, it makes sense to compare this to ROMS-native-reader.

[arian-dialectaquiz]

I did as you suggested, and compared the get_variables () class with the respective ROMS and adjusted obsolete parts in the ECOM reader, rewrote some steps and performed some tests. The main one was with the sigma interpolation for Z, but that also didn't show any change.

The problem, however, remains the same.

The reader_ECOM_S2Z with its reviews is attached below.

From where come the xgrid and ygrid? Directly from the netcdf file? reader_ECOM_S2Z.txt

[gauteh]

In structured there is __make_projected_grid__ and __validate_projected_grid__ which might be useful in checking whether you have valid coordinate arrays.

https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/basereader/structured.py#L437 https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/basereader/structured.py#L469

[knutfrode]

I would also suggest to insert some print-statements in get_variables in your reader. E.g. you may print some variable indices (indx, indy, ...) at various stages to see whether these look meaningful.

For datasets where projection is not known (such as ROMS, and your model), xgrid and ygrid correspond to the row- and column indices (could also be called i and j) of the arrays in the dataset. 2D-polynomial interpolators are then made and used to convert back and forth between (lon, lat) <-> (i, j)

For readers where projection is known, Proj library is used to convert back and forth between (lon, lat) <-> (x, y). Proj uses exact equations, and is faster and more accurate than using polynomial interpolators.

[arian-dialectaquiz]

Thank you! This last comment gave a good clarification. I believe me and my colleague have identified the problem because of it. As soon as I try what we have in mind, I will return here.

But indx and indy remain throughout the stages of the reader, always representing the position in lat / lon -> y / x of the points required over time.

[arian-dialectaquiz]

Goog afternoon! I did as suggested and tried to look at how the projections are in the dataset.

In the ECOM dataset, xgrid and ygrid are xpos and ypos, so I renamed them to be xgrid and ygrid.

And then xgrid and ygrid came into existence as a projection.

Such xgrid and ygrid have minimum and maximum.

However, the problem remains the same. Apparently I still haven't managed to get interpolators.py to receive the xgrid from the dataset. In the Linear2D method the xgrid and ygrid remain the same as indx and indy,respectively. Where indx is the longitude correspondent and indy the latitude correspondent after lat / lon -> x / y interpolation

[knutfrode]

You are already providing 2D arrays of lon and lat, and this is all OpenDrift needs. The screenshot of print(reader) above shows that the geolocation is already perfectly fine. So you should not provide any xgrid/ygrid arrays in addition, these are not used for anything.

The problem is probably instead in your get_variables method, where wrong results are returned, for some reason. I suggest looking at this method step-by-step, to make sure data and indices are reasonable at each stage. e.g., you can insert print(variables) before this dictionary is returned from the reader at the last line. Or you may e.g. plot these data with matplotlib, for debugging.

[arian-dialectaquiz]

Hallo!

I performed step-by-step tests within the variables method, comparing it with the ROMS reader and the basereader's variables.py.

The only differences from reader_ROMS_native are the sigma interpolation for z and the buffer.

The interpolation from sigma to z, even ignored - forcing the reader to add all variables on the surface (indz = 0) - does not change the problem, and does not add any more problem.

But adding a buffer does solve it.

Now xgrid max and min are different:

And the buffer size problem, that was small for the case was solved as well.

Thank you very much, again for your help and willingness!

[knutfrode]

Very good! Then I am closing this issue for now. Please reopen if there are new problems.