#2554 closed help (fixed)

Regriddng standard grid to non-standard squashy grid needed by NEMO

Reported by: charlie Owned by: um_support
Component: UM Model Keywords:
Cc: Platform: NEXCS
UM Version: 10.7

Description

Hi,

Please would you be able to offer some advice on the easiest way to regrid a standard grid (I believe the term is rectilinear?) at 1° resolution to a non-standard variable grid (not sure if the term is curvilinear or unstructured)?

I am trying to regrid a bathymetry file, herold_etal_eocene_topo_1x1.nc (which is global, 1° resolution, standard) to the resolution of the bathymetry file used by NEMO, eORCA_R1_bathy_meter_v2.0.nc (which is global, but has a squashy grid around the tropics). So I already have a source file and a destination file. Both are on JASMIN (although I can easily transfer them to Nexcs, which is where the final version will need to be), at: /group_workspaces/jasmin2/ncas_climate/users/cwilliams2011/data1_for_use.d/sweet.d/gc31.d/ancils.d/bathy.d

I have been told that, because of the squashy grid, I can't use relatively simple tools such as CDO, and instead I need to use either NCL or ESMF. However, I have read the various documentation for both of these but, having never used either, am a bit out of my depth in knowing how to do this and don't really know where to start. Please can you help?

Thanks,

Charlie

Change History (15)

comment:1 Changed 12 months ago by andy

Hi Charlie,

This should be possible in cf-python.
a - herold_etal_eocene_topo_1x1.nc looks like paelo land data with no data over the oceans.
b - eORCA_R1_bathy_meter_v2.0.nc looks like an ocean grid with present day continent outlines.

You'd like to translate the grid from a to that of b. Does that make sense physically as ocean and land are in different places on these grids?

Thanks
Andy

comment:2 Changed 12 months ago by charlie

Hi Andy,

Yes, that's essentially right - I need to take a and translate it to the great of b.

Currently a should have data over both land and ocean, where land has positive values (topography, in metres) and ocean has negative values (bathymetry, again in metres). Given that I am only interested in the bathymetry this time, my plan is to firstly zero everything over land, so I just end up with data over the oceans. This, as you say, uses a paleao mask. I then want to regrid this from his native 1° square resolution to the variable resolution of b. I don't believe it matters that the ocean and land are in different places in the 2 grids, because I don't care about losing the modern mask. All I care about is maintaining the ocean data from the paleao mask, but having it on the variable grid that is required by NEMO. Does that make sense?

Charlie

comment:3 Changed 12 months ago by andy

Hi Charlie,

The following looks to work in cf-python with using ESMF regridding in that it producesa file called cw.c that is an interpolated version of the regular grid onto the orca grid.

python2.7

import cf, cfplot as cfp
regular=cf.read('herold_etal_eocene_topo_1x1.nc')[0]
orca=cf.read('eORCA_R1_bathy_meter_v2.0.nc')[0]
new = regular.regrids(orca, method='bilinear', dst_axes={'X': 'ncdim%x', 'Y': 'ncdim%y'})
cf.write(new, 'cw.nc')

The regrids method in cf-python is documented at:
https://cfpython.bitbucket.io/docs/latest/generated/cf.Field.regrids.html

comment:4 Changed 12 months ago by charlie

Thanks very much Andy, I will try that right now. Does it also copy the attributes, or do I need to do that separately? Although obviously the attributes for the field won't match the original anyway, because they have different names.

Charlie

comment:5 Changed 12 months ago by charlie

Hi Andy,

Sorry, but it's not working quite as I needed to. The code worked when I tried it on the original herold data, but not now once I try it on the same file with the positive values zeroed.

You can see my various steps at /home/users/cwilliams2011/analysis.d/projects.d/sweet.d/gc31n96orca1_eo.d/build.d/ocean.d/bathy.d on JASMIN@CEDA. The actual data are at /group_workspaces/jasmin2/ncas_climate/users/cwilliams2011/data1_for_use.d/sweet.d/gc31.d/ancils.d/bathy.d/

Firstly, I use IDL (make1_zerotopo.pro) to read in the original file (herold_etal_eocene_topo_1x1.nc) and make all positive values (which correspond to topography over land) zero. I then write this out to a new netcdf file (herold_etal_eocene_topo_1x1_notopo.nc). This seems to work fine.

Secondly, I used python (make2_getattributes.py) to get the correct attributes in my new file (this code was given to me by Luke for another purpose, but works fine here as well). Rather than copying the attributes from the original into my newly created file, the code does it the other way round i.e. copies the new data into original. So I end up with the original file (herold_etal_eocene_topo_1x1.nc) being overwritten, such that it contains all the correct attributes AND my new data. This also seems to work.

Lastly, I use your python (which I have now named make3_changegrid.py) to regrid this to the ORCA resolution. This, however doesn't work - I get the following error:

[cwilliams2011@jasmin-sci2 bathy.d]$ python2.7 make3_changegrid.py
Segmentation fault (core dumped)

and, although the brand-new file is created, xconv doesn't recognise it.

Has the code actually worked after all, and xconv doesn't recognise it because it's the wrong version (which was a problem on Monsoon, meaning I had to use a different version of xconv) or has the code not worked?

Thanks,

Charlie

comment:6 Changed 12 months ago by andy

Hi Charlie,

I cannot get into /home/users/cwilliams2011/analysis.d/projects.d/sweet.d/gc31n96orca1_eo.d/build.d/ocean.d/bathy.d to see what's happening here. Could you
cmhod 755 /home/users/cwilliams2011
and then chmod -R 755 /home/users/cwilliams2011/analysis.d

so that I can gain access?

Thanks
Andy

comment:7 Changed 12 months ago by charlie

Sorry Andy, that's caught me several times on Jasmin. Now done.

Charlie

comment:8 Changed 12 months ago by andy

Hi Charlie,

Just tested the final step on jasmin3 - this is the high memory interactive server. This worked fine for me and I can see the file okay on jasmin-sci3 using xconv. I think maybe your conversion ran out of memory on another sci server and the file ended up being corrupted. Can you also do this okay on jasmin-sci3?

Thanks
Andy

comment:9 Changed 12 months ago by charlie

Hi Andy,

Very many thanks, I was using sci2 so clearly that was the problem. If I run it using sci3, it works and create the file no problem.

However, I don't think it's doing quite what it should. If I look at the resulting file, there are several issues: firstly the regridding appears to have gone wrong in certain places (e.g. squares over Antarctica); secondly the map orientation has been changed (such that Africa is now on the right rather than the left); and lastly, and most importantly, it appears to have blanked out the ocean. The land points appear to be near zero (which they should be, or rather exactly zero, because that's exactly what my first step, make1*pro, does), but the ocean points only appear to have data around the edges and not throughout the ocean.

If you compare the original (herold_etal_eocene_topo_1x1_orig.nc) with the newly created file (herold_etal_eocene_topo_orcagrid_notopo.nc), they should be virtually identical over the ocean, with the only difference being the land should be all zero in the new version. As you can see in the original, positive values are over land (and indicate topography) whereas negative values are over ocean (and indicate bathymetry). I want to end up with a file that maintains these negative values, maintains the continents in the same place, but on the ORCA grid and with everything over land zeroed out.

Charlie

comment:10 Changed 12 months ago by andy

Hi Charlie,

This is very confusing and I've spend an hour looking at various files and made no real progress.

I'd like to go back to where we started this discussion and use the files herold_etal_eocene_topo_1x1.nc and eORCA_R1_bathy_meter_v2.0.nc from seven days ago. The file /group_workspaces/jasmin2/ncas_climate/users/cwilliams2011/data1_for_use.d/sweet.d/gc31.d/ancils.d/bathy.d/herold_etal_eocene_topo_1x1.nc has changed and I'd like to get the original. Could you point me to where they are so I can pick the original up?

Thanks
Andy

comment:11 Changed 12 months ago by charlie

Sorry this is causing such a problem Andy, it might be that python can’t do what I need which is why I was advised to use NCL instead.

To answer your question: the original file should be in that same directory but ending _orig.nc however if you want the absolute unadulterated original it is one directory back in herold.d ie cd ../herold.d and then has the same name.

I also realise there is a final step I need to do, once we have the new version of the herold data (ie without positive values over land) and once this has been re gridded to the ORCA resolution, I then need to reinsert this back into the original ORCA file. As you can see from this, there about 7 or 8 fields in this file but all of the others can stay as they are ie it’s only the bathymetry that we are changing. So my final product should look exactly like the original ORCA file, with the same number of fields, but including the new bathymetry that is based on the herold file but having removed positive values over land.

Does that make sense?

Charlie

comment:12 Changed 12 months ago by andy

Hi Charlie,

I think you'll need to go down the NCL route here. There are some anomalies in the Nemo data that are causing artefacts in the cf-python interpolated data. This will need a fix from ESMF in the US and that might take some time to come through.

Cheers
Andy

comment:13 Changed 12 months ago by charlie

Hi Andy,

Okay. In that case, given that I have never used NCL, would you be able to help with using this to do what I need to do? I have been told that the NCL regridding instructions are at https://www.ncl.ucar.edu/Applications/ESMF.shtml but, to be honest, I don’t really understand which script I should be using? I don’t even know how to start NCL or indeed which platform I should be using?!

Charlie

comment:14 Changed 12 months ago by andy

Hi Charlie,

I've never used ncl so I'm not going to be any help here. ncl is available on jasmin@CEDA just by typing ncl.

Cheers
Andy

comment:15 Changed 11 months ago by andy

  • Resolution set to fixed
  • Status changed from new to closed

We went to email to sort this issue out.

There's some interesting things going on here.

1) xconv ignores any masking when making plots
2) cf-python honours masking

Your data herold_etal_eocene_topo_1x1.nc has a _FillValue of -999 but no valid_min, valid_max or valid_range. It does have a actual_range which I've not heard of before.

On https://github.com/Unidata/netcdf4-python/issues/748 there's a discussion of what the netCDF4 Python library does in such a circumstance:

From the netcdf users guide under valid_range:

"If neither valid_min, valid_max nor valid_range is defined then generic applications should define a valid range as follows. If the data type is byte and _FillValue is not explicitly defined, then the valid range should include all possible values. Otherwise, the valid range should exclude the _FillValue (whether defined explicitly or by default) as follows. If the _FillValue is positive then it defines a valid maximum, otherwise it defines a valid minimum. For integer types, there should be a difference of 1 between the _FillValue and this valid minimum or maximum. For floating point types, the difference should be twice the minimum possible (1 in the least significant bit) to allow for rounding error."

So, the code is interpreting the_FillValue as a valid_min. I agree this may be surprising, but that's why I did it.

To get around this use the nco tool ncatted to insert valid_min and valid_max values
cp herold_etal_eocene_topo_1x1.nc herold2.nc
ncatted -a valid_min,topo,c,f,-10000.0 herold2.nc
ncatted -a valid_max,topo,c,f,10000.0 herold2.nc

In cf-python we can now see that the data looks better:

# Old data
f=cf.read('herold_etal_eocene_topo_1x1.nc')[0]
print f.min(), f.max()
# gives -998.6015625 m 2730.16186523 m

# New data
f=cf.read('herold2.nc')[0]
print f.min(), f.max()
# gives -5305.45263672 m 2730.16186523 m

So now we can get started on some cf-python code:

import cf
import numpy as np

# Read in
f=cf.read('herold2.nc')[0]

# Set all values ≥0 to be 0
g = f.where(f≥0, 0)

# Check data looks okay
print g.min(), g.max()
## gives -5305.45263672 m 0.0 m

# Nemo expects bathymetry to be positive values
# so lets multiply by -1
g = g * -1

# Check data looks okay
print g.min(), g.max()
## gives -0.0 m 5305.45263672 m

# regrid to Orca and write out the data
orca=cf.read('eORCA_R1_bathy_meter_v2.0.nc')[0]
new = g.regrids(orca, method='bilinear', dst_axes={'X': 'ncdim%x', 'Y': 'ncdim%y'})
cf.write(new, 'bathymetry_orca.nc')

Note: See TracTickets for help on using tickets.