Adding projection to ROMS native reader

[tklenz]

Hi there,

I am trying to calculate LCS's from a ROMS simulation and am running into a bit of an issue. Using the native reader "as is" the fakeproj correctly returns the corners of my domain, but only places one particle in it.

===========================
Reader: roms native
Projection:
None
Coverage: [pixels]
xmin: 0.000000 xmax: 793.000000 step: 1 numx: 794
ymin: 0.000000 ymax: 361.000000 step: 1 numy: 362
Corners (lon, lat):
(199.67, 55.23) (213.36, 63.39)
(206.19, 52.26) (220.02, 59.70)
[...]

09:27:29 INFO opendrift.models.basemodel:1701: All points are in ocean
09:27:29 INFO opendrift.models.basemodel:2870: 2016-08-02 16:00:00 - step 1 of 120 - 1 active elements (0 deactivated)
09:27:29 INFO opendrift.models.basemodel:2870: 2016-08-02 15:50:00 - step 2 of 120 - 1 active elements (0 deactivated)

Consequently, no ftles can be calculated. From a previous thread (issue #783) I gathered that I need to add a projection to the reader and did so in opendrift.readers.basereader.variables line 14:

class ReaderDomain(Timeable):
"""
Projection, spatial and temporal domain of reader.
"""
name = None

proj4 = '+proj=lcc +lat_1=55 +lat_2=60 +lon_0=-95'
proj = None

Using print(reader) I get this:

===========================
Reader: roms native
Projection:
+proj=lcc +lat_1=55 +lat_2=60 +lon_0=265 +units=m
Coverage: [degrees]
xmin: -3531750.000000 xmax: -2342250.000000 step: 1500 numx: 794
ymin: -3003750.000000 ymax: -2462250.000000 step: 1500 numy: 362
Corners (lon, lat):
(-111.88, -15.26) (-106.32, -14.12)
(-111.27, -17.62) (-105.91, -16.57)

and finally particles are placed on the grid

17:00:50 INFO opendrift.models.basemodel:2762: Backwards simulation, starting at time of last seeded element
17:01:45 INFO opendrift.models.basemodel:1689: Using existing reader for land_binary_mask
17:01:45 INFO opendrift.models.basemodel:1701: All points are in ocean
17:01:45 INFO opendrift.models.basemodel:2870: 2016-08-02 16:00:00 - step 1 of 120 - 644980 active elements (0 deactivated)
17:01:53 INFO opendrift.models.basemodel:2870: 2016-08-02 15:50:00 - step 2 of 120 - 644980 active elements (0 deactivated)

The problem is that the corners of my domain as displayed by print(reader) are not correct (compare to using fakeproj above). How do I fix this? Do I need to provide the projection info somewhere else?

Thank you!

[knutfrode]

Hi,

Inserting projection information inside the reader code will have uncontrolled implication, and here leads to error.

I believe instead the problem is that the delta input parameter to calculate_ftle us using native coordinates, which for ROMS native is not in meters (as in examples found in the gallery), but rather in pixels.
Thus if your model pixel size is e.g. 4 km, you should use delta=1 to have one particle per ROMS pixel, og e.g. delta=3 to place particles with spacing of 12 km.

Note that the LCS/FTLE methods are quite experimental, and should probably be moved outside of OpenDrift, as a wrapper/postprocessor functionality.
But you can anyway try if adjusting delta will give meaningful output.

[simonweppe]

Hi @knutfrode
I have a question about ROMS projection as well (but not related to FTLE).

If used "as is" the reader will use the fakeproj that assumes:

No proj string or projection could be derived, using 'fakeproj'. This assumes that the variables are structured and gridded 
approximately equidistantly on the surface (i.e. in meters). This must be guaranteed by the user. You can get rid of this warning by 
supplying a valid projection to the reader.

In our case ROMS grid is provided as 2D arrays (like the example netcdf file) lon_rho, lat_rho . If I understand well, if the grid is such that grid cell dimensions are consistent across the grid, then it's fine to use that fakeproj, right ?
I havent dived back in the code, but that basically means that velocity at particles positions will be interpolated just linearly using the 2D arrays lon/lat correct ? (i.e. find closest nodes and do a distance weighting for the interpolated value?). Are you using that for the Norway scale runs with ROMS ?

If grid cell dimensions are NOT consistent across the grid, it seems we cant actually provide a projection to the ROMS reader even if we wanted to. Passing proj4 when initializing wont work but I suppose it can be done like this, after initializing the reader e.g.

nordic_native = reader_ROMS_native.Reader(o.test_data_folder() +
    '2Feb2016_Nordic_sigma_3d/Nordic-4km_SLEVELS_avg_00_subset2Feb2016.nc') # proj4 = '+proj=lonlat +ellps=WGS84')
# try to add the projection info afterward ?
import pyproj
nordic_native.proj4 = '+proj=lonlat +ellps=WGS84'
nordic_native.proj = pyproj.Proj(nordic_native.proj4)

Keen to hear your thoughts

[knutfrode]

Hi,
That warning is one of many in OpenDrift that sounds more scary than it actually is :-)
Yes, when projection is not known, bilinear interpolation of 2D arrays of lon and lat is used to convert between (i,j) and (lon,lat).

I believe this is fine for all practical purposes.
Perhaps we should make this warning sound a bit less scary, @gauteh ?

[simonweppe]

Thanks for the clarification @knutfrode

[kthyng]

@knutfrode I'm still confused about whether I should be inputting a proj4 string to the ROMS reader. It looks like if I were to do this, then it would save time calculating the "interpolator for lon,lat to x,y conversion", which I would like to avoid doing to save time, if possible. However, I cannot input "proj4" to the ROMS reader. Should I add this as an optional input to the ROMS reader or was it intentionally removed as an option?

[knutfrode]

According to the docstring it seems one can provide proj4-string to the reader, but this is a recent mistake when the docstring was copied from reader_netCDF_CF_generic:
45251a3

If a proj4-string were to be used, the file would also need to contain X- and Y-coordinates. But I believe in this case it would anyway be better to use reader_netCDF_CF_generic (and provide proj4 there, if the projection if not detected automatically).
I believe the ROMS-native reader would need to be re-written a little if it should use 1D X/Y-coordinates instead of 2D lon/lat arrays.

Note that there is some legacy and historical mess in the ROMS native reader. One day it will hopefully be merged with the generic reader.
It was written first time some years ago when ROMS native output was not very CF-compliant. However, I guess in the meantime this has become better(?). Personally I have not used the ROMS native reader for some years, as we internally produce regridded/post-processed files automatically.

[kthyng]

@knutfrode Thanks for the info. I do have lon/lat arrays as well as x/y arrays for each grid. What do you think is the best way to use these in the ROMS reader so they don't have to be calculated? I understand now that using reader_netCDF_CF_generic for ROMS output is probably preferable in the long run, but probably that is more of a long term goal.

[kthyng]

@knutfrode Alternatively, I thought I could alter the projection code in structured.py to calculated a Delaunay triangulation to input to LinearNDInterpolator which could be reused across simulations to save time.

[kthyng]

@knutfrode Ok it turns out I can pickle the interpolators you already have in there to save time. I made changes to this effect in a new PR (I'm about to submit).

One question though: the "proj" error says that the grid points should be approximately equidistant — is that true? In structured.py it looks like it is accepting the unstructured (as in, unraveled and not equidistant) lon/lat locations to make the interpolator between lon/lat and x/y. However, I want to make sure because my model output is curvilinear and it is using this interpolation approach.

[knutfrode]

Yes, curvilinear should be just fine, grid points are then "approximately equidistant", and there should be only very minor (insignificant) numerical errors in interpolation compared to a regular grid (which anyway only strictly applies to a map projection approximation).

This warning might sound a bit more scary than it is, and should perhaps be removed or "mildened". What do you think @gauteh ?

[kthyng]

@knutfrode Sorry to keep pushing on this, but given that LinearNDInterpolator accepts the lon/lat values from the model and they can have whatever spacing between them, doesn't that mean that the interpolators are accurately representing the conversion between lon/lat and index space for the model grid (at least up to the accuracy of the Qhull representation calculated in LinearNDInterpolator)?