Project

General

Profile

Shifting Longitudes

Added by Markus Ritschel almost 3 years ago

I am having a problem with one of the grids of the GFDL-CM4 model. It seems like the longitudes of the areacello gn grid span a range from -300...60.
Is there a way to shift the longitudes via CDO such that they either range from -180...+180 or from 0...+360?
Thanks a lot in advance!
Markus


Replies (8)

RE: Shifting Longitudes - Added by Markus Ritschel almost 3 years ago

Btw, I also tried `sellonlatbox,-180,180,-90,90` as suggested here: https://code.mpimet.mpg.de/boards/1/topics/22
But without success...

RE: Shifting Longitudes - Added by Karin Meier-Fleischer almost 3 years ago

Hi Markus,

can you upload the file (one timestep)?

-Karin

RE: Shifting Longitudes - Added by Markus Ritschel almost 3 years ago

Hi Karin,

Thanks a lot for your reply. Please find the file appended.

Cheers,
Markus

RE: Shifting Longitudes - Added by Karin Meier-Fleischer almost 3 years ago

As always ;-), the curvilinear grid is the problem, and I guess that you have to interpolate the data to a regular grid.

RE: Shifting Longitudes - Added by Markus Ritschel almost 3 years ago

Those blo***, wonderful curvilinear grids :-D
But there is no way according to your knowledge that retains the original grid? Because I would like to avoid regridding my data as much as possible.

RE: Shifting Longitudes - Added by Karin Meier-Fleischer almost 3 years ago

Hi Markus,

a workaround is to use NCL to read the file content and shift lon and lon_bnds by 360. (=> lon min=60.0005 max=420.034). The lon and lon_bnds have to be shifted like that to make sure that the curvilinear grid won't be corrupted.

NCL script shift_lon.ncl:

f = addfile("GFDL-CM4_areacello_gn.nc","r") 

areacello      = f->areacello
lon2d    = f->lon
lat2d    = f->lat
lat_bnds = f->lat_bnds
lon_bnds = f->lon_bnds
xh       = f->xh
yh       = f->yh

lon2d_shift    = lon2d + 360.
lon_bnds_shift = lon_bnds + 360.

copy_VarAtts(lon2d, lon2d_shift)
copy_VarAtts(lon_bnds, lon_bnds_shift)

;-- write shifted lon2d_shift and lon_bnds_shift to netCDF file

outfile = "$HOME/Downloads/GFDL-CM4_areacello_gn_lon_shift.nc" ;-- output file name

if(isfilepresent(outfile)) then    
   system("rm -rf "+outfile)                        ;-- make sure that outfile does not exist
end if    

fout = addfile(outfile, "c")                        ;-- open file for creation

;-- begin output file settings    
setfileoption(fout, "DefineMode", True)             ;-- explicitly declare file definition mode

;-- create global attributes of the file
fAtt                  =  True                       ;-- assign file attributes
fAtt@title            = "NCL Efficient Approach to netCDF Creation"  
fAtt@source_file      = "GFDL-CM4_areacello_gn.nc" 
fAtt@Conventions      = "CF"   
fAtt@creation_date    =  systemfunc ("date")        
fAtt@history          = "NCL script: shift_lon_curvilinear.ncl" 
fAtt@comment          = "Shift longitude by +360."       
fileattdef(fout, fAtt)                              ;-- copy file attributes    

;-- get the dimension names and sizes for x and y
dimNames = (/"yh", "xh", "y",  "x",  "vertex"/)     ;-- dimension names
dimSizes = (/1080, 1440, 1080, 1440, 4/)            ;-- dimension sizes
dimUnlim = (/ False, False, False, False, False/)   ;-- dimension types
filedimdef(fout, dimNames, dimSizes, dimUnlim)

;-- predefine the the dimensionality of the variables to be written out
filevardef(fout, "xh",        typeof(xh),        getvardims(xh)) 
filevardef(fout, "yh",        typeof(yh),        getvardims(yh)) 
filevardef(fout, "lat",       typeof(lat2d),     getvardims(lat2d)) 
filevardef(fout, "lon",       typeof(lon2d),     getvardims(lon2d)) 
filevardef(fout, "areacello", typeof(areacello), getvardims(areacello)) 
filevardef(fout, "lat_bnds",  typeof(lat_bnds),  getvardims(lat_bnds)) 
filevardef(fout, "lon_bnds",  typeof(lon_bnds),  getvardims(lon_bnds)) 

;-- copy attributes associated with each variable to the file
filevarattdef(fout,"xh",        xh)
filevarattdef(fout,"yh",        yh)
filevarattdef(fout,"lat",       lat2d)
filevarattdef(fout,"lon",       lon2d_shift)
filevarattdef(fout,"lat_bnds",  lat_bnds) 
filevarattdef(fout,"lon_bnds",  lon_bnds_shift)
filevarattdef(fout,"areacello", areacello)

;-- explicitly exit file definition mode (not required)
setfileoption(fout,"DefineMode",False)

;-- Output only the data values since the dimensionality and such have been predefined.
fout->xh        =  (/xh/)
fout->yh        =  (/yh/)
fout->lat       =  (/lat2d/)
fout->lon       =  (/lon2d_shift/)
fout->lat_bnds  =  (/lat_bnds/)
fout->lon_bnds  =  (/lon_bnds_shift/)
fout->areacello =  (/areacello/)

delete(fout)

Run NCL script:

ncl shift_lon.ncl
> ncdump -h GFDL-CM4_areacello_gn_lon_shift.nc

netcdf GFDL-CM4_areacello_gn_lon_shift {
dimensions:
    yh = 1080 ;
    xh = 1440 ;
    y = 1080 ;
    x = 1440 ;
    vertex = 4 ;
variables:
    double xh(xh) ;
        xh:long_name = "h point nominal longitude" ;
        xh:units = "degrees_east" ;
        xh:cartesian_axis = "X" ;
    double yh(yh) ;
        yh:long_name = "h point nominal latitude" ;
        yh:units = "degrees_north" ;
        yh:cartesian_axis = "Y" ;
    float lat(yh, xh) ;
        lat:long_name = "latitude" ;
        lat:units = "degrees_north" ;
        lat:missing_value = 1.e+20f ;
        lat:_FillValue = 1.e+20f ;
        lat:cell_methods = "y: point x: point" ;
        lat:standard_name = "latitude" ;
        lat:bounds = "lat_bnds" ;
        lat:axis = "Y" ;
    float lon(yh, xh) ;
        lon:missing_value = 1.e+20f ;
        lon:axis = "X" ;
        lon:bounds = "lon_bnds" ;
        lon:standard_name = "longitude" ;
        lon:cell_methods = "y: point x: point" ;
        lon:units = "degrees_east" ;
        lon:long_name = "longitude" ;
        lon:_FillValue = 1.e+20f ;
    float areacello(yh, xh) ;
        areacello:_FillValue = 1.e+20f ;
        areacello:cell_methods = "area: sum" ;
        areacello:long_name = "Grid-Cell Area" ;
        areacello:missing_value = 1.e+20f ;
        areacello:standard_name = "cell_area" ;
        areacello:units = "m2" ;
        areacello:coordinates = "lat lon" ;
        areacello:original_name = "areacello" ;
    float lat_bnds(y, x, vertex) ;
        lat_bnds:long_name = "latitude bounds" ;
        lat_bnds:units = "degrees_north" ;
        lat_bnds:axis = "Y" ;
    float lon_bnds(y, x, vertex) ;
        lon_bnds:long_name = "longitude bounds" ;
        lon_bnds:units = "degrees_east" ;
        lon_bnds:axis = "X" ;

// global attributes:
        :comment = "Shift longitude by +360." ;
        :history = "NCL script: shift_lon_curvilinear.ncl" ;
        :creation_date = "Mi  5 Mai 2021 15:22:23 CEST" ;
        :Conventions = "CF" ;
        :source_file = "GFDL-CM4_areacello_gn.nc" ;
        :title = "NCL Efficient Approach to netCDF Creation" ;
}

The output file is attached.

-Karin

RE: Shifting Longitudes - Added by Markus Ritschel almost 3 years ago

Hi Karin,

Wow, thanks a lot for this extensive answer! I will check it out later. Thanks also for appending the processed file.

All the best,
Markus

RE: Shifting Longitudes - Added by Karin Meier-Fleischer almost 3 years ago

At the moment the script above will only shift lon and lon_bnds, not xh. Maybe you need to do the same with xh, it depends on the application you want to use further on.

    (1-8/8)