target_grid_cube for iris-esmf-regrid

[jmpmcmanus]

What is the best way to create a target_grid_cube file for iris-esmf-regrid? I had a geotiff that I converted to netcdf, but it thru and error because it did not have bounds. I added bounds using ncap2:

ncap2 -s ' defdim("bnds",2); lon_bnds=make_bounds(lon,$bnds); lat_bnds=make_bounds(lat,$bnds);' maxwvel.raw.63.nc maxwvel.test.63.nc

but it thru and error on that file as well.

Thanks
Jim

[stephenworsley]

Hi Jim,

Have you tried loading the cube into iris and calling cube.coord("latitude").guess_bounds() and cube.coord("longitude").guess_bounds()? What do your cubes and coordinates look like when printed before regridding?

[jmpmcmanus]

Hi Stephen,

That worked, but when I try to save the data to a netcdf file:

iris.save(result_cube, '/Users/jmpmcman/Work/Surge/data/apsvizvolume/regrid/maxwvel.regrid.63.nc')

I get the following error message:

ValueError: Expected an array whose shape ends in (3564104,), got an array with shape ending in (1813443,).

This same error occurs when I print(type(result.data)).

The shape of the result_cube is (8000, 8838), which is the same as the target_grid_cube, which I was expecting The values of 3564104, and 1813443, are from the source_mesh_cube, which has:

node = 1813443 ;
nele = 3564104 ;

I'm not sure why those values appear in the regridded result_cube? I don't see them in any of the result_cube meta data.

The output from an ncdump -c on the netCDF file used for the target_grid_cube, minus the actual lat/lon values is:

netcdf maxwvel.raw.63 {
dimensions:
lon = 8838 ;
lat = 8000 ;
variables:
char crs ;
crs:grid_mapping_name = "latitude_longitude" ;
crs:long_name = "CRS definition" ;
crs:longitude_of_prime_meridian = 0. ;
crs:semi_major_axis = 6378137. ;
crs:inverse_flattening = 298.257223563 ;
crs:spatial_ref = "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]" ;
crs:crs_wkt = "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]" ;
crs:GeoTransform = "-98.0021427904 0.005000049052613714 0 45.838806 0 -0.005000111564074999 " ;
double lat(lat) ;
lat:standard_name = "latitude" ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
double lon(lon) ;
lon:standard_name = "longitude" ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
double Band1(lat, lon) ;
Band1:long_name = "GDAL Band Number 1" ;
Band1:_FillValue = NaN ;
Band1:grid_mapping = "crs" ;

// global attributes:
:GDAL_AREA_OR_POINT = "Area" ;
:Conventions = "CF-1.5" ;
:GDAL = "GDAL 3.4.2, released 2022/03/08" ;
:history = "Wed May 18 13:50:14 2022: GDAL CreateCopy( maxwvel.raw.63.nc, ... )" ;
data:

and ncdump -c on the netCDF files used for the source_mesh_cube is:

netcdf maxwvel.63 {
dimensions:
time = UNLIMITED ; // (1 currently)
node = 1813443 ;
nele = 3564104 ;
nvertex = 3 ;
nope = 1 ;
neta = 186 ;
max_nvdll = 186 ;
nbou = 186 ;
nvel = 62972 ;
max_nvell = 28573 ;
mesh = 1 ;
variables:
double time(time) ;
time:long_name = "model time" ;
time:standard_name = "time" ;
time:units = "seconds since 2020-12-15 00:00:00" ;
time:base_date = "2020-12-15 00:00:00" ;
double x(node) ;
x:long_name = "longitude" ;
x:standard_name = "longitude" ;
x:units = "degrees_east" ;
x:positive = "east" ;
double y(node) ;
y:long_name = "latitude" ;
y:standard_name = "latitude" ;
y:units = "degrees_north" ;
y:positive = "north" ;
int element(nele, nvertex) ;
element:long_name = "element" ;
element:cf_role = "face_node_connectivity" ;
element:start_index = 1 ;
element:units = "nondimensional" ;
int adcirc_mesh(mesh) ;
adcirc_mesh:long_name = "mesh_topology" ;
adcirc_mesh:cf_role = "mesh_topology" ;
adcirc_mesh:topology_dimension = 2 ;
adcirc_mesh:node_coordinates = "x y" ;
adcirc_mesh:face_node_connectivity = "element" ;
int neta ;
neta:long_name = "total number of elevation specified boundary nodes" ;
neta:units = "nondimensional" ;
int nvdll(nope) ;
nvdll:long_name = "number of nodes in each elevation specified boundary segment" ;
nvdll:units = "nondimensional" ;
int max_nvdll ;
int ibtypee(nope) ;
ibtypee:long_name = "elevation boundary type" ;
ibtypee:units = "nondimensional" ;
int nbdv(neta) ;
nbdv:long_name = "node numbers on each elevation specified boundary segment" ;
nbdv:units = "nondimensional" ;
int nvel ;
nvel:long_name = "total number of normal flow specified boundary nodes including both the front and back nodes on internal barrier boundaries" ;
nvel:units = "nondimensional" ;
int nvell(nbou) ;
nvell:long_name = "number of nodes in each normal flow specified boundary segment" ;
nvell:units = "nondimensional" ;
int max_nvell ;
int ibtype(nbou) ;
ibtype:long_name = "type of normal flow (discharge) boundary" ;
ibtype:units = "nondimensional" ;
int nbvv(nvel) ;
nbvv:long_name = "node numbers on normal flow boundary segment" ;
nbvv:units = "nondimensional" ;
double depth(node) ;
depth:long_name = "distance below geoid" ;
depth:standard_name = "depth below geoid" ;
depth:coordinates = "time y x" ;
depth:location = "node" ;
depth:mesh = "adcirc_mesh" ;
depth:units = "m" ;
double wind_max(node) ;
wind_max:long_name = "maximum wind velocity" ;
wind_max:standard_name = "maximum_wind" ;
wind_max:units = "m s-1" ;
wind_max:_FillValue = -99999. ;
wind_max:coordinates = "y x" ;
wind_max:location = "node" ;
wind_max:mesh = "adcirc_mesh" ;
double time_of_wind_max(node) ;
time_of_wind_max:long_name = "time of maximum wind velocity" ;
time_of_wind_max:standard_name = "time_of_maximum_wind" ;
time_of_wind_max:units = "sec" ;
time_of_wind_max:_FillValue = -99999. ;
time_of_wind_max:coordinates = "y x" ;
time_of_wind_max:location = "node" ;
time_of_wind_max:mesh = "adcirc_mesh" ;

// global attributes:
:_FillValue = -99999. ;
:model = "ADCIRC" ;
:version = "v53.04-17-ga629de8-modified" ;
:grid_type = "Triangular" ;
:description = "HSOFS cs:20201215000000 cy:2022042506 ASGS NAM ! 32 CHARACTER ALPHANUMERIC RUN" ;
:agrid = "NOMAD mesh v1e MSL" ;
:rundes = "HSOFS cs:20201215000000 cy:2022042506 ASGS NAM ! 32 CHARACTER ALPHANUMERIC RUN" ;
:runid = "84 hour namforecast run ! 24 CHARACTER ALPANUMERIC RUN IDENTIFICATION" ;
:title = "ASGS HSOFS V1.0.1 ! NCPROJ - PROJECT TITLE" ;
:institution = "Seahorse Coastal Consulting, LLC ! NCINST - PROJECT INSTITUTION" ;
:source = "ADCIRC Surge Guidance System (ASGS) ! NCSOUR - PROJECT SOURCE" ;
:history = "CSDL/CMMB - NWS/NHC ! NCHIST - PROJECT HISTORY" ;
:references = "http://www.nauticalcharts.noaa.gov/csdl/RD.html ! NCREF - PROJECT REFERENCES" ;
:comments = "HSOFS V1.0.1 ! NCCOM - PROJECT COMMENTS" ;
:host = "https://github.com/jasonfleming/asgs ! NCHOST - PROJECT HOST" ;
:convention = "CF-1.1 ! NCCONV - CONVENTIONS" ;
:Conventions = "UGRID-0.9.0" ;
:contact = "[email protected] ! NCCONT - CONTACT INFORMATION" ;
:creation_date = "2022-04-25 4:44:58 -04:00" ;
:modification_date = "2022-04-25 4:44:58 -04:00" ;
:fort.15 = "==== Input File Parameters (below) ====" ;
:dt = 2. ;
:ihot = 568 ;
:ics = 2 ;
:nolibf = 1 ;
:nolifa = 2 ;
:nolica = 0 ;
:nolicat = 0 ;
:nwp = 5 ;
:ncor = 1 ;
:ntip = 1 ;
:nws = -12 ;
:nramp = 1 ;
:tau0 = -3. ;
:statim = 0. ;
:reftim = 0. ;
:rnday = 499.75 ;
:dramp = 5. ;
:a00 = 0. ;
:b00 = 1. ;
:c00 = 0. ;
:h0 = 0.05 ;
:slam0 = -71. ;
:sfea0 = 27. ;
:cf = 0.0025 ;
:eslm = 10. ;
:cori = 0. ;
:ntif = 8 ;
:nbfr = 8 ;
data:

time = 43178400 ;
}

Thanks!
Jim

[stephenworsley]

By the looks of things, I think the problem here is that the source_mesh_cube is defined on the node location rather than the face location. Area weighted conservative regridding requires data to be located on faces. With that said, the next version of iris-esmf-regrid ought to have bilinear regridding which allows data to be located on nodes. This should also come with better error handling so that this kind of problem is easier to detect.

[jmpmcmanus]

That's good to hear! When will the next version be out? Do you know of any short term work arounds?

Saving the regridder also throughs an error:

File /opt/anaconda3/envs/pangeo/lib/python3.8/site-packages/iris/experimental/ugrid/mesh.py:2283, in _Mesh1DCoordinateManager.filter(self, **kwargs)
2279 name = "" if name is None else f"{name!r} "
2280 emsg = (
2281 f"Expected to find exactly 1 {name}coordinate, but found none."
2282 )
-> 2283 raise CoordinateNotFoundError(emsg)
2285 return result

CoordinateNotFoundError: 'Expected to find exactly 1 coordinate, but found none.'

Thanks!
Jim

[stephenworsley]

I would imagine that the regridder not saving properly is due to it trying to create a MeshCoord on the face location. This attempts to create a Coord with bounds but no points which iris doesn't currently allow. This is a something we will be looking to fix, as per #4438 this is likely to be a priority for UGRID development. In the mean time, it should be possible to work around this by artificially putting in face coordinates, though I think that if you were trying to save a bilinear regridder I think the saving process would be creating MeshCoords on node locations so you would avoid this issue.

I can't give you a date for when the next version of iris-esmf-regrid will be out, though we're trying to keep updates fairly frequent so I'll see if we can get a release out soon enough.

[jmpmcmanus]

Hi Stephen - Thanks again for you help on this! I will look into defining the mesh on the face location, instead of the node location. Do you know if Iris has an option like guess_bounds (target_grid_cube.coord("latitude").guess_bounds()), that you showed me, for defining the mesh on the face location?

[stephenworsley]

I don't believe there's any convenience functions we have to do this yet. It looks like you'd have to go through a couple steps to manually manipulate the data. First of all, you would need to add latitude and longitude face coords to the mesh via mesh.add_coords(face_x=..., face_y=...). This should probably be an average of the coordinates of the nodes, but for area weighted regridding this can just be a zero valued dummy coordinate and regridding ought to still work. Then you will need a new cube defined on the coordinates mesh.to_MeshCoords("face"). bear in mind this will have a different shape to your original cube and you will have to do some kind of interpolation onto this which I don't think we have a shortcut for either. I suppose another approach would be to replace the mesh with its dual (i.e. replace each face center with a node and node with a face center) which is also likely to be very complicated, though wouldn't require interpolating the data. So I'd say this is possible but quite awkward to do currently.