opendrift icon indicating copy to clipboard operation
opendrift copied to clipboard

Sandard name mapping when using the unstructured grid reader

Open IGanch opened this issue 3 years ago • 11 comments

Hello, I am trying to run a simulation using the unstructured grid reader. I want to use the standard_name_mapping because the reader doesn't recognize the names notation, but I get an error for unexpected keyword argument. Is there a way to work around this problem? I am changing the names with nco, but so far the changed netCDF file doesn't look as expected. Many thanks in advance

IGanch avatar Apr 08 '22 08:04 IGanch

Hi,

The standard_name_mapping is presently only supported by the generic and ROMS readers, but should be straightforward to also implement for the unstructured reader. I will have not have time to implement and test this whithin the coming weeks, but in principle you only need to copy (possibly adapt) a few lines from the genric reader into reader_netCDF_CF_unstructured.py: https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_netCDF_CF_generic.py#L117 https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_netCDF_CF_generic.py#L327 If you make it to work, a git pull request is welcome.

If this works, you should not need to manipulate your netCDF files with nco (which is an alternative).

knutfrode avatar Apr 08 '22 09:04 knutfrode

Thanks a lot for the quick reply. I believe my assumption that the unstructured reader would solve my problem was wrong. I am using forcing data for the currents from NEMO and the standard NetCDF reader doesn't recognize the area that the data is covering and gives wrong corners for the data coverage. Can the corners of the data be assigned as well?

IGanch avatar Apr 08 '22 10:04 IGanch

No, the geolocation has to be obtained from the file/URL. What is the error message that you get? Maybe you can post the ncdump of this file?

knutfrode avatar Apr 08 '22 10:04 knutfrode

The error message is ValueError: Readers must be added for the following required variables: ['x_sea_water_velocity', 'y_sea_water_velocity'] Here's the NetCDFs info:

netcdf ASIS1_1h_20180707_20180707_grid_U {
dimensions:
        y = 261 ;
        x = 591 ;
        nvertex = 4 ;
        depthu = 121 ;
        axis_nbounds = 2 ;
        time_counter = UNLIMITED ; // (24 currently)
variables:
        float area(y, x) ;
                area:standard_name = "cell_area" ;
                area:units = "m2" ;
        float bounds_nav_lat(y, x, nvertex) ;
        float bounds_nav_lon(y, x, nvertex) ;
        float depthu(depthu) ;
                depthu:name = "depthu" ;
                depthu:long_name = "Vertical U levels" ;
                depthu:units = "m" ;
                depthu:positive = "down" ;
                depthu:bounds = "depthu_bounds" ;
        float depthu_bounds(depthu, axis_nbounds) ;
                depthu_bounds:units = "m" ;
        float nav_lat(y, x) ;
                nav_lat:standard_name = "latitude" ;
                nav_lat:long_name = "Latitude" ;
                nav_lat:units = "degrees_north" ;
                nav_lat:bounds = "bounds_nav_lat" ;
        float nav_lon(y, x) ;
                nav_lon:standard_name = "longitude" ;
                nav_lon:long_name = "Longitude" ;
                nav_lon:units = "degrees_east" ;
                nav_lon:bounds = "bounds_nav_lon" ;
        double time_centered(time_counter) ;
                time_centered:standard_name = "time" ;
                time_centered:long_name = "Time axis" ;
                time_centered:calendar = "gregorian" ;
                time_centered:units = "seconds since 1900-01-01 00:00:00" ;
                time_centered:time_origin = "1900-01-01 00:00:00" ;
                time_centered:bounds = "time_centered_bounds" ;
        double time_centered_bounds(time_counter, axis_nbounds) ;
        double time_counter(time_counter) ;
                time_counter:axis = "T" ;
                time_counter:standard_name = "time" ;
                time_counter:long_name = "Time axis" ;
                time_counter:calendar = "gregorian" ;
                time_counter:units = "seconds since 1900-01-01 00:00:00" ;
                time_counter:time_origin = "1900-01-01 00:00:00" ;
                time_counter:bounds = "time_counter_bounds" ;
        double time_counter_bounds(time_counter, axis_nbounds) ;
        float uo(time_counter, depthu, y, x) ;
                uo:standard_name = "sea_water_x_velocity" ;
                uo:long_name = "ocean current along i-axis" ;
                uo:units = "m/s" ;
                uo:online_operation = "average" ;
                uo:interval_operation = "200 s" ;
                uo:interval_write = "1 h" ;
                uo:cell_methods = "time: mean (interval: 200 s)" ;
                uo:cell_measures = "area: area" ;
                uo:_FillValue = 1.e+20f ;
                uo:missing_value = 1.e+20f ;
                uo:coordinates = "time_centered nav_lat nav_lon" ;

// global attributes:
                :name = "ASIS1_1h_20180707_20180707_grid_U" ;
                :description = "ocean U grid variables" ;
                :title = "ocean U grid variables" ;
                :Conventions = "CF-1.6" ;
                :timeStamp = "2021-Oct-10 20:51:33 GMT" ;
                :uuid = "02cd57ca-c7e1-4b32-b782-1adebc92acb6" ;
                :history = "Wed Oct 20 09:14:57 2021: ncks -O -7 -L 5 a20180717/ASIS1_1h_20180707_20180707_grid_U.nc a20180717/ASIS1_1h_20180707_20180707_grid_U.nc" ;
                :NCO = "netCDF Operators version 4.8.1 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
}

The x and y current velocities are in tow separate file, but I read in a previous post that this should not be a problem. Many thanks.

IGanch avatar Apr 08 '22 10:04 IGanch

This looks ok, and this file should be read with reader_netCDF_CF_generic What do you get from this?

>>> from opendrift.readers import reader_netCDF_CF_generic
>>> r = reader_netCDF_CF_generic.Reader('ASIS1_1h_20180707_20180707_grid_U')
>>> print(r)

knutfrode avatar Apr 08 '22 10:04 knutfrode

The output is

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 suppling a valid projection to the reader.
===========================
Reader: ASIS1_1h_20180707_20180707_grid_U.nc
Projection:
  None
Coverage: [pixels]
  xmin: 0.000000   xmax: 590.000000   step: 1   numx: 591
  ymin: 0.000000   ymax: 260.000000   step: 1   numy: 261
  Corners (lon, lat):
    ( -1.00,  -1.00)  ( -1.00,  -1.00)
    ( 27.26,  40.50)  ( -1.00,  -1.00)
Vertical levels [m]:
  Not specified
Available time range:
  start: 2018-07-07 00:30:00   end: 2018-07-07 23:30:00   step: 1:00:00
    24 times (0 missing)
Variables:
  cell_area
  latitude
  longitude
  x_sea_water_velocity
===========================

IGanch avatar Apr 08 '22 10:04 IGanch

This also looks ok, except that some of the corner coordinates are clearly wrong (as you commented). It seems that the lon- lat arrays (nav_lon, nav_lat) are corrupted, or as I remember from an earlier issue, perhaps they are not defined over land(?) Thus I believe the only fix for this would be to make the lon- and lat arrays complete. You may plot them e.g. with ncview to see how they look.

knutfrode avatar Apr 08 '22 10:04 knutfrode

The lat lon are ok as I inspected them. The issue is that they cover an area, which is not rectangular and on the empty spaces the values are filled with -1. Here's the plot for longitude. This is the Black Sea basin. image I thought a workaround would be to make a rectangular subset which covers a field without empty values, but this turned out also not so straight forward so far. Many thanks.

IGanch avatar Apr 08 '22 10:04 IGanch

Yes, the problem is that OpenDrift expects the lon- and lat arrays to be complete 2D-arrays without holes like here. OpenDrift cannot understand that the -1 are not real positions, and thus some corners get the coordinates (-1, -1)

I see that @ocgabs had the same problem in this issue: https://github.com/OpenDrift/opendrift/issues/550 Perhaps ocgabs found a workaround in the meantime?

knutfrode avatar Apr 08 '22 11:04 knutfrode

Ok. So a successful subset would be a solution. Perhaps ocgabs has another idea as well. Thanks a lot for your work.

IGanch avatar Apr 08 '22 11:04 IGanch

Yes, if you can make a subset (e.g. with ncks) with only ocean pixels, it should work.

If this is a standard output from NEMO (since you seem both have this kind of file independently of each other?), we could consider later to make a fix in the netCDF generic reader for this special case. Although we would like to avoid too many such hacks in OpenDrift.

knutfrode avatar Apr 08 '22 11:04 knutfrode