Project

General

Profile

Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates

Added by S BR about 3 years ago

Hi All,
I have a regional gridded data in the British National Grid (OSGB) spatial coordinate system. I am finding to use this data in CDO? Can we at least regrid this data to Latitude-longitude in rotated pole coordinates or to a regular lat/lon grid, so that I can do other analysis using CDO.

For your convenience I have attached here a sample NetCDF file.

Thank you.
SBR


Replies (13)

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Brendan DeTracey about 3 years ago

See section 1.4.2.4 of the documentation. Also, a search for British National Grid gave this topic:
https://code.mpimet.mpg.de/boards/1/topics/10198
Does that answer your question?

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by S BR about 3 years ago

I have tried as per the given information but there is still no success. I couldn't generate the correct file.

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Brendan DeTracey about 3 years ago

Without your command line (with -v verbose option), your grid description file, and the output text, it is almost impossible to help you. I understand why you can't upload your very large data file, but could you at least ncdump -h your_netcdf_file.nc and post the output text? Not having any information on your data file makes help doubly impossible. :)

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by S BR about 3 years ago

Hi Brendan,
I have already uploaded the file in my first post. Again I uploading here. At first I checked with cdo sinfo tas_rcp85_land-rcm_uk_12km.nc and you can see that the file shows wrong coordinate information (I have informed that this file plots correctly in Python as this file is generated with Iris python)
File format : NetCDF4 classic
    -1 : Institut Source   T Steptype Levels Num    Points Num Dtype : Parameter ID
     1 : unknown  UKCP18   v instant       1   1      9184   1  F32  : -1
     2 : unknown  UKCP18   v instant       1   2         1   2  I32  : -2
   Grid coordinates :
     1 : curvilinear              : points=9184 (82x112)
                   grid_longitude : -18.24146 to -6.580251 degrees_east
                    grid_latitude : -0.8710619 to 12.94984 degrees_north
                          mapping : transverse_mercator
          projection_x_coordinate : -210000 to 762000 by 12000 m
          projection_y_coordinate : -102000 to 1230000 by 12000 m
                        available : cellbounds
     2 : generic                  : points=1
   Vertical coordinates :
     1 : generic                  : levels=1
                  ensemble_member : 1 1
     2 : surface                  : levels=1
   Time coordinate :  10 steps

If I do ncdump -h, I get the following information.
netcdf tas_rcp85_land-rcm_uk_12km {
dimensions:
    time = UNLIMITED ; // (10 currently)
    bnds = 2 ;
    projection_x_coordinate = 82 ;
    projection_y_coordinate = 112 ;
    ensemble_member = 1 ;
variables:
    double time(time) ;
        time:standard_name = "time" ;
        time:bounds = "time_bnds" ;
        time:units = "hours since 1970-01-01 00:00:00" ;
        time:calendar = "360_day" ;
        time:axis = "T" ;
    double time_bnds(time, bnds) ;
    double grid_longitude(projection_y_coordinate, projection_x_coordinate) ;
        grid_longitude:standard_name = "longitude" ;
        grid_longitude:long_name = "longitude" ;
        grid_longitude:units = "degrees_east" ;
        grid_longitude:_CoordinateAxisType = "Lon" ;
    double grid_latitude(projection_y_coordinate, projection_x_coordinate) ;
        grid_latitude:standard_name = "latitude" ;
        grid_latitude:long_name = "latitude" ;
        grid_latitude:units = "degrees_north" ;
        grid_latitude:_CoordinateAxisType = "Lat" ;
    double projection_x_coordinate(projection_x_coordinate) ;
        projection_x_coordinate:standard_name = "projection_x_coordinate" ;
        projection_x_coordinate:units = "m" ;
        projection_x_coordinate:axis = "X" ;
        projection_x_coordinate:bounds = "projection_x_coordinate_bnds" ;
    double projection_x_coordinate_bnds(projection_x_coordinate, bnds) ;
    double projection_y_coordinate(projection_y_coordinate) ;
        projection_y_coordinate:standard_name = "projection_y_coordinate" ;
        projection_y_coordinate:units = "m" ;
        projection_y_coordinate:axis = "Y" ;
        projection_y_coordinate:bounds = "projection_y_coordinate_bnds" ;
    double projection_y_coordinate_bnds(projection_y_coordinate, bnds) ;
    int transverse_mercator ;
        transverse_mercator:grid_mapping_name = "transverse_mercator" ;
        transverse_mercator:longitude_of_prime_meridian = 0. ;
        transverse_mercator:semi_major_axis = 6377563.396 ;
        transverse_mercator:semi_minor_axis = 6356256.909 ;
        transverse_mercator:longitude_of_central_meridian = -2. ;
        transverse_mercator:latitude_of_projection_origin = 49. ;
        transverse_mercator:false_easting = 400000. ;
        transverse_mercator:false_northing = -100000. ;
        transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ;
    double ensemble_member(ensemble_member) ;
        ensemble_member:long_name = "ensemble_member" ;
        ensemble_member:units = "1" ;
        ensemble_member:axis = "Z" ;
    float tas(time, ensemble_member, projection_y_coordinate, projection_x_coordinate) ;
        tas:standard_name = "air_temperature" ;
        tas:long_name = "Mean air temperature" ;
        tas:units = "degC" ;
        tas:grid_mapping = "transverse_mercator" ;
        tas:coordinates = "grid_latitude grid_longitude" ;
        tas:_FillValue = 1.e+20f ;
        tas:missing_value = 1.e+20f ;
        tas:description = "Mean air temperature" ;
        tas:label_units = "°C" ;
        tas:plot_label = "Mean air temperature at 1.5m (°C)" ;
        tas:cell_methods = "time: mean" ;
    int year(time) ;
        year:long_name = "year" ;
        year:units = "1" ;

// global attributes:

As per your suggestions w.r.t the post https://code.mpimet.mpg.de/boards/1/topics/10198, I tried the following way. I don't know where to get the the 'proj_params' information, so copied directly from the above post. The I did the following way.
cdo remapnn,grid.txt tas_rcp85_land-rcm_uk_12km.nc newdata.nc
where my grid.txt is given below. Still the new file couldn't produce the correct lat/lon grid.

gridtype = projection
xsize = 82
ysize = 112
xunits   = "meter" 
yunits   = "meter" 
xfirst = -210000
xinc = 12000
yfirst = -102000
yinc = 12000
grid_mapping_name = crs
proj_params = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs" 

Thanks.

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Brendan DeTracey about 3 years ago

Sorry mate. My first response was from misreading your question. My second response was meant for another thread. I'll blame lack of sleep.
cdo should allow you to reproject this to a curvilinear lon-lat grid, but even after a couple hours of fussing I can't do it. Maybe one of the cdo gods can help.

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Ralf Mueller about 3 years ago

hi guys!

I checked the interpolation to a lonlat grid, but the area is shifted from UK to east Africa. So I think CDO has problems with this horizontal grid. A little cross-check seems to uncover, that the internal grid representation is not self-sufficient in this case:

  1. get rid of the year variable (non-CF-conform):
    cdo selname,tasmax tasmax_rcp85_land-rcm_uk_12km_time10.nc tasmax_proj.nc
  2. write out the internal grid description
    cdo griddes tasmax_proj.nc > gd
  3. perform nearest neighbour interpolation onto itself
    cdo -L remapnn,gd tasmax_proj.nc again_tasmax_proj.nc

    the last step throws an error:
    cdo    remapnn (Abort): nvertex undefined (grid description file: gd)!

So I think in the current state, CDO cannot correctly work with this input grid, sorry.

best wishes
ralf

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Karin Meier-Fleischer about 3 years ago

Hi,

and now, I weigh in here, too. ;)

We also had the discussion at the ncl-talk mailing list, and the two data files were not remapable and still stays over Africa. But there is a solution. :)

To remap the data to a lonlat grid over GB, first create a grid file (target grid) using the topo operator:

cdo -f nc -sellonlatbox,-14,6,48,62 -remapbil,r1440x720 -topo topo.nc

Get the grid information of the input file using the griddes operator. You have to modify the grid information file (see uploaded grid_OSGB.txt file).
@Brendan DeTracey: You were close!

Then remap your input file to the lonlat target grid:

cdo -remapbil,topo.nc -setgrid,grid_OSGB.txt -selname,tas tas_rcp85_land-rcm_uk_12km.nc outfile.nc

It is important to use the setgrid oprator, too.

-Karin

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Brendan DeTracey about 3 years ago

I will add my last two pfennigs. Since you know how to handle your grid in python, I suggest you could create a cdo grid file description(https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-210001.4.2) (https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-910000D.1) in longitude and latitude. Corner points can be easily created in projection_x_coordinate/projection_y_coordinate space with constant offsets from the center points and then reprojected using pyproj or equivalent. It will take a little time and effort but will give you optimal results. setgrid your data to this curvilinear grid and the end result is that you do not need to interpolate at all. Your data file coordinates will be real lat and lon on the original grid. It is always better to try to do all your processing on the original grid. After that if you need to remap your processed results, you can easily do so.

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Karin Meier-Fleischer about 3 years ago

@Brendan DeTracey: Your're right, doing calculations on the origin grid is the best, but the issue was to remap the input data to a lonlat grid. Most applications, also the python packages, do the grid handling under the hood which means you do not have the original latitude/longitude values of the grid. The grid_OSGB.txt grid description file now contains the needed grid information.

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Brendan DeTracey about 3 years ago

@Karin. I guess I am over-complicating again. My idea was that if someone did have knowledge of python/pyproj/numpy, they could easily create a curvilinear grid in projection_coordinates in meters, transform these coordinates to true latitude and longitude, and write to a grid description, if the original file was CF compliant(Section 5.6). Not as easy, but possible.

Sounds like this thread would be a good addition to the Operator News topic. Don't look at me! ;)

RE: Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates - Added by Karin Meier-Fleischer about 3 years ago

@Brendan DeTracey: Yes, of course that would make a good Operator News article. If I ever have some time... ;)

    (1-13/13)