Project

General

Profile

Forums » Support »

Projection parameter is missing for remapping

Added by Simon Werner 12 months ago

Dear Support-community:
I have this file with climate data

"netcdf dSM_map {
dimensions:
time = UNLIMITED ; // (2557 currently)
easting = 49 ;
northing = 42 ;
variables:
int crs ;
crs:proj4 = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs" ;
int time(time) ;
time:standard_name = "time" ;
time:long_name = "time" ;
time:units = "hours since 1951-01-01 00:00:00" ;
time:calendar = "standard" ;
time:axis = "T" ;
double lon(northing, easting) ;
lon:standard_name = "longitude" ;
lon:long_name = "longitude" ;
lon:units = "deg. dec." ;
lon:_CoordinateAxisType = "Lon" ;
double lat(northing, easting) ;
lat:standard_name = "latitude" ;
lat:long_name = "latitude" ;
lat:units = "deg. dec." ;
lat:_CoordinateAxisType = "Lat" ;
double easting(easting) ;
easting:standard_name = "projection_x_coordinate" ;
easting:long_name = "x-coordinate in the given coordinate system" ;
easting:units = "m or deg. dec." ;
easting:axis = "X" ;
double northing(northing) ;
northing:standard_name = "projection_y_coordinate" ;
northing:long_name = "y-coordinate in the given coordinate system" ;
northing:units = "m or deg. dec." ;
northing:axis = "Y" ;
float SWC ;
SWC:coordinates = "lat lon" ;
SWC:_FillValue = -9999.f ;
SWC:missing_value = -9999.f ;
SWC:grid_mapping = "crs" ;
SWC:proj4 = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs" ;

// global attributes:
:CDI = "Climate Data Interface version 2.0.5 (https://mpimet.mpg.de/cdi)" ;
:Conventions = "CF-1.6" ;
:title = "mHMv5.5 simulation outputs" ;
:creation_date = "2016-06-21 13:20:20" ;
:CDO = "Climate Data Operators version 2.0.5 (https://mpimet.mpg.de/cdo)" ;
:NCO = "netCDF Operators version 5.0.5 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
:history = "Fri Jan 5 14:23:23 2024: ncatted -a proj4,SWC,o,c,+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs tmp.nc dSM_map.nc\nFri Jan 5 14:23:21 2024: ncatted -a grid_mapping,SWC,o,c,crs dSM_thu_time_crs.nc tmp.nc\nFri Jan 5 13:42:34 2024: ncatted -O -a proj4,crs,o,c,+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs temp.nc dSM_thu_time_crs.nc\nFri Jan 5 13:42:34 2024: ncap2 -O -s crs=0 dSM_thu_time.nc temp.nc\nThu Nov 23 17:51:18 2023: cdo selyear,2016/2022 dSM_thu.nc dSM_thu_time.nc\nThu Nov 23 17:40:53 2023: cdo sellonlatbox,9.88,12.67,50.19,51.66 dSM.nc dSM_thu.nc\nWed Oct 25 15:46:26 2023: cdo mergetime SM_p1.nc SM_p2.nc dSM.nc\nWed Oct 25 15:44:41 2023: cdo -expr,SWC=SWC_L01+SWC_L02 SM1.nc SM_p1.nc\nWed Oct 25 15:17:15 2023: cdo selvar,SWC_L01,SWC_L02 /data/duerre/GDM_v1/data/MHM_r2229/output/mHM_Fluxes_States.nc SM1.nc" ;
}"

And this grid file functioning as my source grid:

"netcdf grid_r_crs {
dimensions:
easting = 50 ;
northing = 41 ;
variables:
int crs ;
crs:proj4 = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs" ;
double easting(easting) ;
easting:units = "m o. deg dec." ;
easting:long_name = "x-coordinate in the given coordinate system" ;
easting:standard_name = "projection_x_coordinate." ;
easting:axis = "X" ;
double northing(northing) ;
northing:units = "m o. deg dec." ;
northing:long_name = "y-coordinate in the given coordinate system" ;
northing:standard_name = "projection_y_coordinate." ;
northing:axis = "Y" ;
float layer(northing, easting) ;
layer:_FillValue = -3.4e+38f ;
layer:grid_mapping = "crs" ;
layer:proj4 = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs" ;
layer:min = 17065. ;
layer:max = 24067. ;

// global attributes:
:Conventions = "CF-1.4" ;
:created_by = "R, packages ncdf4 and raster (version 3.6-23)" ;
:date = "2024-01-05 11:55:43" ;
:history = "Fri Jan 5 14:34:56 2024: ncatted -O -a axis,northing,o,c,Y tmp.nc grid_r_crs.nc\n",
"Fri Jan 5 14:34:56 2024: ncatted -O -a standard_name,northing,o,c,projection_y_coordinate. tmp.nc tmp.nc\n",
"Fri Jan 5 14:34:56 2024: ncatted -O -a units,northing,o,c,m o. deg dec. tmp.nc tmp.nc\n",
"Fri Jan 5 14:34:56 2024: ncatted -O -a long_name,northing,o,c,y-coordinate in the given coordinate system tmp.nc tmp.nc\n",
"Fri Jan 5 14:34:56 2024: ncatted -O -a axis,easting,o,c,X tmp.nc tmp.nc\n",
"Fri Jan 5 14:34:56 2024: ncatted -O -a standard_name,easting,o,c,projection_x_coordinate. tmp.nc tmp.nc\n",
"Fri Jan 5 14:34:56 2024: ncatted -O -a units,easting,o,c,m o. deg dec. tmp.nc tmp.nc\n",
"Fri Jan 5 14:34:55 2024: ncatted -O -a long_name,easting,o,c,x-coordinate in the given coordinate system grid_r.nc tmp.nc" ;
:NCO = "netCDF Operators version 5.0.5 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
}"

I am not able to remap the file although the grid description files are consistent in my view:
"$ cdo griddes dSM_map.nc #
  1. gridID 1 #
    gridtype = curvilinear
    gridsize = 2058
    xsize = 49
    ysize = 42
    xname = lon
    xdimname = easting
    xlongname = "longitude"
    xunits = "deg. dec."
    yname = lat
    ydimname = northing
    ylongname = "latitude"
    yunits = "deg. dec."
    xvals = 9.86126865839591 ...
    yvals = 51.6501848027876 51.6512229277301 ...
    scanningMode = 64
    gridtype = projection
    gridsize = 2058
    xsize = 49
    ysize = 42
    xname = easting
    xlongname = "x-coordinate in the given coordinate system"
    xunits = "m or deg. dec."
    yname = northing
    ylongname = "y-coordinate in the given coordinate system"
    yunits = "m or deg. dec."
    xfirst = 4352000
    xinc = 4000
    yfirst = 5726000
    yinc = -4000
    scanningMode = 64
    grid_mapping = crs
    proj4 = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs"
$ cdo griddes grid_r_crs.nc #
  1. gridID 1 #
    gridtype = projection
    gridsize = 2050
    xsize = 50
    ysize = 41
    xname = easting
    xlongname = "x-coordinate in the given coordinate system"
    xunits = "m o. deg dec."
    yname = northing
    ylongname = "y-coordinate in the given coordinate system"
    yunits = "m o. deg dec."
    xfirst = 4348579.08935821
    xinc = 4000
    yfirst = 5724362.27541661
    yinc = -4000
    scanningMode = 64
    grid_mapping = crs
    proj4 = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs"
    cdo griddes: Processed 1 variable [0.00s 13MB].
    "

However my remapping code returns this error and I don't know which projection parameter is missing:

$ cdo remapbil,grid_r_crs.nc dSM_map.nc dSM_44gr.nc

cdo remapbil (Abort): Projection parameter missing!


Replies (2)

RE: Projection parameter is missing for remapping - Added by Estanislao Gavilan 12 months ago

Hi Simon,

try to write proj_params instead of proj4. Something like this

proj_params = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +units=m +no_defs"

Lets see how it goes.

RE: Projection parameter is missing for remapping - Added by Simon Werner 12 months ago

This worked for me. Thanks allot!

    (1-2/2)