Menu

Regridding netCDF

Help
2013-12-04
2020-08-26
  • Maria Catala

    Maria Catala - 2013-12-04

    Hello all,

    I have a netcdf file with one variable called Z500 (geopotential height at 500 hPa) and three dimensions (latitude, longitude and time). The longitude size is 26 (from 28 to 60 by 2), the latitude size is 17 (from -20 to 30 by 2) and the time size is 205996 (150 years with 6 hourly measurements).

    My grid has a spatial resolution of 2 degrees and I want to regrid it in order to have a spatial resolution of 2.5 degrees. The new latitude should have a size of 13 (from 30 to 60 by 2.5) and the longitude should have a size of 19 (from -17.5 to 27.5 by 2.5).

    In order to do that I am using ncap2 with the folowing script (called Interpol.nco):

    // Bilinear interpolation
    // Included by rgr.sh

    defdim("latn",13);
    defdim("lonn",19);
    latn[$latn] = {30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60};
    lonn[$lonn] = {-17.5, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20 ,22.5, 25, 27.5};

    out[$time,$latn,$lonn]=0.0;

    // Bi-linear interpolation
    Z500f = bilinear_interp_wrap(Z500,out,latn,lonn,lat,lon);

    After that, in the console:

    ncap2 -0 -S Interpol.nco Hgt_500_1871-2011_3d.nc test.nc

    This is the error message that I get:

    ncap2: ERROR bilinear_interp_wrap(): Dimension size mismatch with input variables

    I thought it was because of the time dimension so I tried to write the following code into the script:
    i=1
    until ($i == dim(time))
    do Z500f[$i,,] = bilinear_interp_wrap(Z500[$i,,],out[$i,,],latn,lonn,lat,lon);
    $i = $i + 1
    done

    But I don't know if it is a good idea or the good syntax.

    Can someone help me?

    Thanks!

     

    Last edit: Maria Catala 2013-12-04
  • henry Butowsky

    henry Butowsky - 2013-12-04

    Hi Maria,
    you are nearly there: try the following script.

    /******* interpol.nco **********/
    defdim("latn",13);
    defdim("lonn",19);
    latn[$latn] = {30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60};
    lonn[$lonn] = {-17.5, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20 ,22.5, 25, 27.5};

    out[$latn,$lonn]=0.0;
    Z500_new[$time,$latn,$lonn]=0.0;

    for( idx=0 ; idx<$time.size ; idx++){
    Z500_new(idx,:,:)=bilinear_interp(Z500(idx,:,:),out,latn,lonn,lat,lon);
    }

    /*******************/

    Then run the script as follows
    ncap2 -v -S interpol.nco in.nc out.nc

    ... Henry

     
    • Craig Cobb

      Craig Cobb - 2016-03-16

      Hi Henry, was looking at example and trying to figure out if this would help. Am trying to take a netcdf file and subset AND regrid to another interval. Can I do this with one command (using NCO version 3 operators)? Or can I build the resampled gridfile I need and use it in a comman? Have looked through examples and haven't found anything. Thanks in advance.

       
  • Maria Catala

    Maria Catala - 2013-12-05

    Thanks for your answer! This is exactly what I tried to do (but I am only used to R language)

    I got several problems, this is what I have with ncdump -h in.nc:

    dimensions:

    lon = 26 ;
    lat = 17 ;
    time = UNLIMITED ; // (5844 currently)
    

    variables:

    float lon(lon) ;
        lon:standard_name = "longitude" ;
        lon:long_name = "Longitude" ;
        lon:units = "degrees_east" ;
        lon:axis = "X" ;
    
    float lat(lat) ;
        lat:standard_name = "latitude" ;
        lat:long_name = "Latitude" ;
        lat:units = "degrees_north" ;
        lat:axis = "Y" ;
    
    double time(time) ;
        time:standard_name = "time" ;
        time:units = "hours since 1800-01-01 00:00:00" ;
    
        time:calendar = "standard" ;
    float Z500(time, lat, lon) ;
        Z500:standard_name = "geopotential_height" ;
        Z500:long_name = "6-hourly Geopotential Heights on Pressure Levels" ;
        Z500:units = "m" ;
        Z500:_FillValue = -32767.f ;
        Z500:unpacked_valid_range = -1500.f, 60000.f ;
        Z500:actual_range = 4866.f, 5998.f ;
        Z500:precision = 0s ;
        Z500:GRIB_id = 7s ;
        Z500:GRIB_name = "HGT" ;
        Z500:var_desc = "Geopotential height" ;
        Z500:dataset = "NOAA-CIRES 20th Century Reanalysis version 2" ;
        Z500:level_desc = "Pressure Levels" ;
        Z500:statistic = "Ensemble Mean" ;
        Z500:parent_stat = "Other" ;
    

    This is what I have with the outpout file:

    dimensions:

    latn = 13 ;
    lonn = 19 ;
    time = UNLIMITED ; // (5844 currently)
    lat = 17 ;
    lon = 26 ;
    

    variables:

    double Z500_new(time, latn, lonn) ;
    
    double latn(latn) ;
    
    double lonn(lonn) ;
    
    double out(latn, lonn) ;
    
    int idx ;
    
    float lon(lon) ;
    
        lon:standard_name = "longitude" ;
        lon:long_name = "Longitude" ;
        lon:units = "degrees_east" ;
        lon:axis = "X" ;
    
    float lat(lat) ;
        lat:standard_name = "latitude" ;
        lat:long_name = "Latitude" ;
        lat:units = "degrees_north" ;
        lat:axis = "Y" ;
    
    double time(time) ;
        time:standard_name = "time" ;
        time:units = "hours since 1800-01-01 00:00:00" ;
        time:calendar = "standard" ;
    

    In fact I want to find the same format as the input file. lonn, latn and Z500 are considered as double and not as float, can I change that ? Moreover Z500_new[,,1] is a grid with 0 (after this time step I have normal values).

     

    Last edit: Maria Catala 2013-12-05
  • henry Butowsky

    henry Butowsky - 2013-12-05

    Hi Maria,
    your editing your previous posts - its getting confusing. Can just reply/post instead of editing your posts.

    The existence of the tmp file indicates that the script has run incorrectly.
    Can you try doing a bilineariterp on the first time step in the series. Just to check things ate working :

    /***********/
    defdim("latn",13);
    defdim("lonn",19);
    // change to your new range here
    latn[$latn] = {30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60};
    // change to your new range here
    lonn[$lonn] = {-17.5, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20 ,22.5, 25, 27.5};
    out[$latn,$lonn]=0.0;

    // no time dimension
    Z500_new[$latn,$lonn]=0.0;

    Z500_new=bilinear_interp(Z500(0,:,:),out,latn,lonn,lat,lon);

    /***/

    ... Henry

     

    Last edit: henry Butowsky 2013-12-05
  • Maria Catala

    Maria Catala - 2013-12-05

    I am sorry Henry, I didn't see your post!

    I had 150 years data so it really took time. I used another file to test (with only three years data) and it works well. The weird thing is when I add the line "print(Z500_new)" the "for" never ends ! But when the line is not here, it works well...

    The output I get is the following:

    dimensions:

    latn = 13 ;
    
    lonn = 19 ;
    
    time = UNLIMITED ; // (5844 currently)
    
    lat = 17 ;
    
    lon = 26 ;
    

    variables:

    double Z500_new(time, latn, lonn) ;
    
    double latn(latn) ;
    
    double lonn(lonn) ;
    
    double out(latn, lonn) ;
    
    int idx ;
    
    float lon(lon) ;
        lon:standard_name = "longitude" ;
        lon:long_name = "Longitude" ;
        lon:units = "degrees_east" ;
        lon:axis = "X" ;
    
    float lat(lat) ;
        lat:standard_name = "latitude" ;
        lat:long_name = "Latitude" ;
        lat:units = "degrees_north" ;
        lat:axis = "Y" ;
    
    double time(time) ;
        time:standard_name = "time" ;
        time:units = "hours since 1800-01-01 00:00:00" ;
        time:calendar = "standard" ;
    

    I used ncrcat -v var -x out.nc out2.nc to delete lon, lat, idx and out variables
    I also used ncwa -a dim out2.nc out3.nc to delete lon and lat dimensions. After that I have:

    dimensions:

    latn = 13 ;
    
    lonn = 19 ;
    
    time = UNLIMITED ; // (5844 currently)
    

    variables:

    double Z500_new(time, latn, lonn) ;
    
    double latn(latn) ;
    
    double lonn(lonn) ;
    
    double time(time) ;
        time:standard_name = "time" ;
        time:units = "hours since 1800-01-01 00:00:00" ;
        time:calendar = "standard" ;
    

    I will use this file in a program and usually it is float or short variables. But I think it will be ok with a double type.

    Now I will try with the "big" file, I'll tell in one hour...

     

    Last edit: Maria Catala 2013-12-05
    • Maria Catala

      Maria Catala - 2013-12-05

      Wow in fact, when I do it with 3 years it takes 30 seconds. But when I do it with 158 years, I've calculated that it will take 17 hours!! I don't know why there is such a difference... I'll let the function run this night...

      Edit: It works well thanks but is there a way to keep the type of the variables? (Short or float for coordinates)

       

      Last edit: Maria Catala 2013-12-06
  • henry Butowsky

    henry Butowsky - 2016-03-17

    HI Craig,
    you can subset and regrid at the same time -
    it may take a little more time if you dont regrid as for each new gridpoint the bounds in the coordinate variables have to be searched
    ....Henry

     
  • henry Butowsky

    henry Butowsky - 2016-03-17

    hi Craig,
    sorry I missunderstood your question.
    With ncap2 you cannot subset a region like you can in the other operators using command line parameters. However if have already prepared the new coordinate vars then there will be no need t osubset
    ....Henry

     
  • Craig Cobb

    Craig Cobb - 2016-03-17

    Hi Henry, I'm sorry but I don't understand your last statement about already having prepared the new coordinate vars... Could you elaborate, particularly with which functions to use. Thanks

    Craig

     
  • henry Butowsky

    henry Butowsky - 2016-03-18

    hi Craig,
    coordinate variables describe a mathematical space 2D/3D /4D
    Typically we use latitude/longitude to describe the surface of the earth -
    Lat /lon on comes in different flavours regular/irregular 1D/2D

    So for example NOAA GSLnetcdf files comes
    with lat - -90.0, -89.5,-89.0...... 89.5, 90.5
    lon 0.0,0.5,1.0,2.0...359.5

    We may want to re-grid to a higher resolution or lower resolution or change the wrapping.
    To create regular new coordinate variables can use the array function in ncap2.

    e.g

    /*****************************************************************/
    defdim("lat_new", 21);
    defdim("lon_new",40);
    lat_new=array(-90.0, 9.0, $lat_new);
    lon_new=array(0.0, 9.0,$lon_new);
    
    // create new temerature variable
    Temp_new[lat_new,lon_new]=0.0f;
    
    // bilinear interpolation
    Temp_new=bilinear_interp(Temp_old,Temp_new,lat_new,lon_new,lat,lon);
    
    /*****************************************************************/
    

    ...Henry

     
    • Craig Cobb

      Craig Cobb - 2016-03-21

      Hi Henry. I was able to download the version 4.5.5 nco operators. Can I regrid my (Netcdf version 3 containing multiple parameters) data to subset and oversample the input with any single command using the new software? Am not familiar with the SCRIP format specified...

      I tried using ncap2 and ncremap with a regridded (map file) but having no luck. I might not have everything set correctly however in the remapped file.

      I'd like to keep it as simple as possible. Thanks.

      Sorry for the questions but am not very familiar with these operators.

       

      Last edit: Craig Cobb 2016-03-21
  • henry Butowsky

    henry Butowsky - 2016-03-23

    Hi Craig,
    to subset use ncks e.g
    ncks -d latitude,-45.0,30.0 -d longitude, 0.0,40.0 in.nc out.nc
    to regrid
    use the bilinear_interpolation as mentioned above
    ....Henry

     
  • Craig Cobb

    Craig Cobb - 2016-03-23

    Thanks. I understand how to subset with ncks, but when using the bilinear_interpolation method doesn't it then target and output single (defined) parameters to another file? Was trying to perform operation over netcdf file containing multiple parameters. Would I add definitions of other variables and perform similar loops? There are 3 and 4 dimensional parameters.

    I used:

    /*/

    defdim("latn",10);
    defdim("lonn",10);
    // change to new latitude range here (using subset of actual values):
    latn[$latn] = {30.046682,30.046960,30.047238,30.047515,30.047793,30.048071,30.048349,30.048626,30.048904,30.049182};
    // change to new longitude range here (using subset of actual values):
    lonn[$lonn] = {274.173284,274.173561,274.173839,274.174116,274.174394,274.174671,274.174948,274.175226,274.175503,274.175781};
    out[$latn,$lonn]=0.0;
    surf_el_new[$time,$latn,$lonn]=0.0;
    for( idx=0 ; idx<$time.size ; idx++ ) {
    surf_el_new(idx,:,:)=bilinear_interp(surf_el(0,:,:),out,latn,lonn,lat,lon);
    }

    /**/

    When I run using:

    ncap2 -v -S interpol.nco /models/RNCOM/amseas_u/netcdf/ncom_relo_amseas_u_2016032300_t024.nc out.nc

    I get the following error:

    ncap2: ERROR assign(): LHS cast for latn - cannot make RHS ~zz@value_list conform.

    Any suggestions to make it work?

    My input netcdf surf_el info is:

    ncom_relo_amseas_u_2016032300_t024.nc: Title: Water Surface Elevation
    ncom_relo_amseas_u_2016032300_t024.nc: Command: not applicable
    ncom_relo_amseas_u_2016032300_t024.nc: Remark:
    ncom_relo_amseas_u_2016032300_t024.nc: Gridline node registration used
    ncom_relo_amseas_u_2016032300_t024.nc: Grid file format: ns (# 16) GMT netCDF format (short) (COARDS-compliant)
    ncom_relo_amseas_u_2016032300_t024.nc: x_min: 262 x_max: 305.096 x_inc: 0.03333 name: [degrees_east] nx: 1294
    ncom_relo_amseas_u_2016032300_t024.nc: y_min: 5 y_max: 32.0973 y_inc: 0.03333 name: [degrees_north] ny: 814
    ncom_relo_amseas_u_2016032300_t024.nc: z_min: -0.774 z_max: 1.86 name: [meter]
    ncom_relo_amseas_u_2016032300_t024.nc: scale_factor: 0.001 add_offset: 0
    ncom_relo_amseas_u_2016032300_t024.nc: 361495 nodes set to NaN

    The header looks like:

    dimensions:
    lon = 1294 ;
    lat = 814 ;
    depth = 40 ;
    time = UNLIMITED ; // (1 currently)
    variables:
    double lon(lon) ;
    lon:units = "degrees_east" ;
    lon:long_name = "Longitude" ;
    lon:NAVO_code = 2 ;
    lon:valid_range = 262., 305.095703125 ;
    double lat(lat) ;
    lat:units = "degrees_north" ;
    lat:long_name = "Latitude" ;
    lat:NAVO_code = 1 ;
    lat:valid_range = 5., 32.0972900390625;
    double depth(depth) ;
    depth:units = "meter" ;
    depth:long_name = "Depth" ;
    depth:NAVO_code = 5 ;
    depth:positive = "down" ;
    double time(time) ;
    time:units = "hour since 2000-01-01 00:00 UTC" ;
    time:long_name = "Valid Time" ;
    time:NAVO_code = 13 ;
    time:time_origin = "2016-03-23 00:00:00" ;
    short surf_el(time, lat, lon) ;
    surf_el:units = "meter" ;
    surf_el:NAVO_code = 32 ;
    surf_el:_FillValue = -30000s ;
    surf_el:missing_value = -30000s ;
    surf_el:long_name = "Water Surface Elevation" ;
    surf_el:add_offset = 0.f ;
    surf_el:scale_factor = 0.001f ;
    surf_el:positive = "up" ;

     

    Last edit: Craig Cobb 2016-03-23
  • henry Butowsky

    henry Butowsky - 2016-03-31

    Hi Craig,
    sorry for the delay in replying - I dont understand why you are getting this wanrning about lonn as it appears to have been initialized with the correct number of floats.
    Can you post your netCDF file to a common location so that i can try doing the interpolation
    ...Henry

     
  • Dawit

    Dawit - 2017-04-22

    Hi All,

    I face similar challange: ncap2: ERROR assign(): LHS cast for latn - cannot make RHS ~zz@value_list conform

    I was trying to regrid prcipitation to higher resolution. The run the script below using ncap2 command as follows:
    ncap2 -v -S pre_rgrid.nco pre_cru_1901_2014.nc pre_cru_new_1901_2014.nc

    /*****/
    defdim("latn",7);
    defdim("lonn",17);
    // change to your new range here
    latn[$latn] = {30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60};
    // change to your new range here
    lonn[$lonn] = {-17.5, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20 ,22.5, 25, 27.5};
    out[$latn,$lonn]=0.0;

    // no time dimension
    pre_cru_new[$latn,$lonn]=0.0;
    for( idx=0 ; idx<$time.size ; idx++ ) {
    pre_cru_new(idx,:,:)=bilinear_interp(pre_cru_new(0,:,:),out,latn,lonn,lat,lon);
    }

    I would appreciate if you share your experince how this challange is addressed.

    Dawit

     
  • Jayant Pendharkar

    Hi All,
    Me too facing the same error: ncap2: ERROR assign(): LHS cast for latitude - cannot make RHS ~zz@value_list conform.
    when regridding to higher resolution. Is there any workaround?

     
    • Charlie Zender

      Charlie Zender - 2020-08-26

      Henry, is it possible to make bilinear_interp() produce more helpful error messages?

       

Log in to post a comment.