CDO

Hello,

I have tried to perform a conservative remapping procedure of RADAR data (1km resolution) onto a curvilinear model grid (approx. 3km resolution) but it seems like the calculated weights are wrong. I think this is because of the (probably) unknown/unsupported grid type of the RADAR observations.
Here are the following steps I have performed (and some additional information):

 cdo griddes CCLM-0275_grid_bnds.nc > CCLM-0275_griddes.txt    # get grid description file of target grid
 cdo gencon,CCLM-0275_griddes.txt RADKLIM_nomiss.nc gencon_wgt.nc # calculate remapping weights
 cdo remap,CCLM-0275_griddes.txt,gencon_wgt.nc RADKLIM_nomiss.nc RADKLIM_remapcon.nc # perform conservative remapping

The result looks like this (comparing original - 1km res with 3km output res):

During the generation of remapping weights many errors occur including map weights above 1 and below 0 (see image).

cdo gencon: Map weight < 0! grid1idx=982738 grid2idx=140606 nlink=1752516 wts=-3.32649e+15
cdo gencon: Map weight > 1! grid1idx=982739 grid2idx=140192 nlink=1752517 wts=8.96177e+14
[...]
cdo gencon: Error: sum of wts for map1 42616 0.199219 1
[...]

Here are some outputs of the verifygrid and ncdump -h of the input and output grids:

>cdo -verifygrid RADKLIM_nomiss.nc
cdo verifygrid: Grid consists of 990000 cells (type: curvilinear), of which
cdo verifygrid:    990000 cells have 4 vertices
cdo verifygrid:     61996 cells are not unique
cdo verifygrid:    990000 cells have their vertices arranged in a clockwise order
cdo verifygrid: Processed 2 variables ( 0.94s 182MB )

>cdo -verifygrid CCLM-0275_grid_bnds.nc
cdo verifygrid: Grid consists of 180525 cells (type: curvilinear), of which
cdo verifygrid:    180525 cells have 4 vertices
cdo verifygrid: Processed 1 variable ( 0.19s 45MB )

>ncdump -h RADKLIM_nomiss.nc
netcdf RADKLIM_nomiss {
dimensions:
        time = UNLIMITED ; // (1 currently)
        x = 900 ;
        y = 1100 ;
        vertices = 4 ;
variables:
        double time(time) ;
                time:standard_name = "time" ;
                time:long_name = "Time" ;
                time:units = "hours since 2001-01-01 00:50:00.0" ;
                time:calendar = "standard" ;
                time:axis = "T" ;
        double lon(y, x) ;
                lon:standard_name = "longitude" ;
                lon:long_name = "Longitude" ;
                lon:units = "degrees_east" ;
                lon:_CoordinateAxisType = "Lon" ;
                lon:bounds = "lon_bnds" ;
        double lon_bnds(y, x, vertices) ;
        double lat(y, x) ;
                lat:standard_name = "latitude" ;
                lat:long_name = "Latitude" ;
                lat:units = "degrees_north" ;
                lat:_CoordinateAxisType = "Lat" ;
                lat:bounds = "lat_bnds" ;
        double lat_bnds(y, x, vertices) ;
        double x(x) ;
                x:standard_name = "projection_x_coordinate" ;
                x:long_name = "RADOLAN Grid x coordinate of projection" ;
                x:units = "m" ;
                x:axis = "X" ;
        double y(y) ;
                y:standard_name = "projection_y_coordinate" ;
                y:long_name = "RADOLAN Grid y coordinate of projection" ;
                y:units = "m" ;
                y:axis = "Y" ;
        int crs ;
                crs:long_name = "RADOLAN Grid" ;
                crs:grid_mapping_name = "polar_stereographic" ;
                crs:latitude_of_projection_origin = 90. ;
                crs:straight_vertical_longitude_from_pole = 10. ;
                crs:semi_major_axis = 6370040. ;
                crs:semi_minor_axis = 6370040. ;
                crs:earth_radius = 6370040. ;
                crs:standard_parallel = 60. ;
                crs:false_easting = 0. ;
                crs:false_northing = 0. ;
                crs:inverse_flattening = 0. ;
                crs:scale_factor_at_projection_origin = 0.9330127019 ;
                crs:units = "m" ;
                crs:proj4 = "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 +k_0=0.933012701892 +x_0=0 +y_0=0 +a=6370040 +b=6370040 +units=km +no_defs" ;
        float RR(time, y, x) ;
                RR:standard_name = "rainfall_amount" ;
                RR:long_name = "hourly rainfall " ;
                RR:units = "kg m-2" ;
                RR:grid_mapping = "crs" ;
                RR:coordinates = "lat lon" ;
                RR:_FillValue = 999.f ;
                RR:missing_value = 999.f ;

>ncdump -h CCLM-0275_grid_bnds
netcdf CCLM-0275_grid_bnds {
dimensions:
        rlon = 415 ;
        rlat = 435 ;
        nv = 4 ;
variables:
        double rlon(rlon) ;
                rlon:long_name = "longitude in rotated pole grid" ;
                rlon:standard_name = "grid_longitude" ;
                rlon:units = "degrees" ;
                rlon:axis = "X" ;
        double rlat(rlat) ;
                rlat:long_name = "latitude in rotated pole grid" ;
                rlat:standard_name = "grid_latitude" ;
                rlat:units = "degrees" ;
                rlat:axis = "Y" ;
        double lon(rlat, rlon) ;
                lon:long_name = "longitude" ;
                lon:standard_name = "longitude" ;
                lon:units = "degrees_east" ;
                lon:bounds = "lon_vertices" ;
        double lat(rlat, rlon) ;
                lat:long_name = "latitude" ;
                lat:standard_name = "latitude" ;
                lat:units = "degrees_north" ;
                lat:bounds = "lat_vertices" ;
        double lon_vertices(rlat, rlon, nv) ;
                lon_vertices:units = "degrees_east" ;
        double lat_vertices(rlat, rlon, nv) ;
                lat_vertices:units = "degrees_north" ;
        char rotated_pole ;
                rotated_pole:long_name = "coordinates of the rotated North Pole" ;
                rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
                rotated_pole:grid_north_pole_latitude = 39.25 ;
                rotated_pole:grid_north_pole_longitude = -162. ;
        byte dummy(rlat, rlon) ;
                dummy:long_name = "dummy_variable" ;
                dummy:grid_mapping = "rotated_pole" ;
                dummy:coordinates = "lat lon" ;

Does anyone have a suggestion on how to solve this issue?
-- Harald