Project generic grid (km) to lat lon grid

Hi there,

I need some ideas to project or convert a .nc file with generic grid coordinates in km to a lat/lon one.

The file is available here; https://zenodo.org/records/4289959/files/smb_rec.2015-2099.BN_RACMO2.3p2-CESM2-SSP5-85.1km.YY.nc.gz?download=1

Below my cod sinfo:

   File format : NetCDF
    -1 : Institut Source   T Steptype Levels Num    Points Num Dtype : Parameter ID
     1 : unknown  unknown  v instant       1   1   4039200   1  F32  : -1            
   Grid coordinates :
     1 : generic                  : points=4039200 (1496x2700)
                                x : -638956 to 856044 by 1000 km
                                y : -3354596 to -655596 by 1000 km
   Vertical coordinates :
     1 : surface                  : levels=1
   Time coordinate :
                             time : 85 steps
     RefTime =  2015-01-01 00:00:00  Units = days  Calendar = standard  Bounds = true
  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss
  2015-06-30 18:00:00  2016-06-30 18:00:00  2017-06-30 18:00:00  2018-06-30 18:00:00
  2019-06-30 18:00:00  2020-06-30 18:00:00  2021-07-02 00:00:00  2022-06-30 18:00:00
  2023-06-30 18:00:00  2024-06-30 18:00:00  2025-06-30 18:00:00  2026-06-30 18:00:00
  2027-06-30 18:00:00  2028-06-30 18:00:00  2029-06-30 18:00:00  2030-06-30 18:00:00
  2031-06-30 18:00:00  2032-06-30 18:00:00  2033-06-30 18:00:00  2034-06-30 18:00:00
  2035-06-30 18:00:00  2036-06-30 18:00:00  2037-06-30 18:00:00  2038-06-30 18:00:00
  2039-06-30 18:00:00  2040-06-30 18:00:00  2041-06-30 18:00:00  2042-06-30 18:00:00
  2043-06-30 18:00:00  2044-06-30 18:00:00  2045-06-30 18:00:00  2046-06-30 18:00:00
  2047-06-30 18:00:00  2048-06-30 18:00:00  2049-06-30 18:00:00  2050-06-30 18:00:00
  2051-06-30 18:00:00  2052-06-30 18:00:00  2053-06-30 18:00:00  2054-06-30 18:00:00
  2055-06-30 18:00:00  2056-06-30 18:00:00  2057-06-30 18:00:00  2058-06-30 18:00:00
  2059-06-30 18:00:00  2060-06-30 18:00:00  2061-06-30 18:00:00  2062-06-30 18:00:00
  2063-06-30 18:00:00  2064-06-30 18:00:00  2065-06-30 18:00:00  2066-06-30 18:00:00
  2067-06-30 18:00:00  2068-06-30 18:00:00  2069-06-30 18:00:00  2070-06-30 18:00:00
  2071-06-30 18:00:00  2072-06-30 18:00:00  2073-06-30 18:00:00  2074-06-30 18:00:00
  2075-06-30 18:00:00  2076-06-30 18:00:00  2077-06-30 18:00:00  2078-06-30 18:00:00
  2079-06-30 18:00:00  2080-06-30 18:00:00  2081-06-30 18:00:00  2082-06-30 18:00:00
  2083-06-30 18:00:00  2084-06-30 18:00:00  2085-06-30 18:00:00  2086-06-30 18:00:00
  2087-06-30 18:00:00  2088-06-30 18:00:00  2089-06-30 18:00:00  2090-06-30 18:00:00
  2091-06-30 18:00:00  2092-06-30 18:00:00  2093-06-30 18:00:00  2094-06-30 18:00:00
  2095-06-30 18:00:00  2096-06-30 18:00:00  2097-06-30 18:00:00  2098-06-30 18:00:00
  2099-06-30 18:00:00
cdo    sinfo: Processed 1 variable over 85 timesteps [0.02s 20MB]

and my ncdump:

netcdf smb_rec.2015-2099.BN_RACMO2.3p2-CESM2-SSP5-85.1km.YY {
dimensions:
    time = UNLIMITED ; // (85 currently)
    y = 2700 ;
    x = 1496 ;
    bnds = 2 ;
variables:
    float smb_rec(time, y, x) ;
        smb_rec:standard_name = "SMB_reconstructed_from_downscaled_components" ;
        smb_rec:long_name = "SMB reconstructed from downscaled components" ;
        smb_rec:units = "mm w.e. per month" ;
        smb_rec:_FillValue = -1.e+30f ;
        smb_rec:missing_value = -1.e+30f ;
        smb_rec:actual_range = -40.12159f, 260.1555f ;
    float time(time) ;
        time:standard_name = "time" ;
        time:long_name = "time" ;
        time:bounds = "time_bnds" ;
        time:units = "DAYS since 2015-01-01 00:00:00" ;
        time:calendar = "standard" ;
        time:axis = "T" ;
    double time_bnds(time, bnds) ;
    float x(x) ;
        x:long_name = "x" ;
        x:units = "km" ;
        x:axis = "X" ;
    float y(y) ;
        y:long_name = "y" ;
        y:units = "km" ;
        y:axis = "Y" ;

// global attributes:
        :CDI = "Climate Data Interface version 1.9.8 (https://mpimet.mpg.de/cdi)" ;
        :Conventions = "CF-1.6" ;
        :institution = "IMAU (Brice Noel)" ;
        :title = "Daily Surface mass balance field (RACMO2.3)" ;
        :grid = "Map Projection:Polar Stereographic Ellipsoid - Map Reference Latitude: 90.0 - Map Reference Longitude: -39.0 - Map Second Reference Latitude: 71.0 - Map Eccentricity: 0.081819190843 ;wgs84 - Map Equatorial Radius: 6378137.0 ;wgs84 meters - Grid Map Origin Column: 160 - Grid Map Origin Row: -120 - Grid Map Units per Cell: 5000 - Grid Width: 301 - Grid Height: 561" ;
        :netcdf = "4.4.1.1 of Jul  4 2018 17:49:55 $" ;
        :frequency = "mon" ;
        :NCO = "netCDF Operators version 4.8.0 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;
        :CDO = "Climate Data Operators version 1.9.8 (https://mpimet.mpg.de/cdo)" ;
}

As you can see coordinates are "project" in km. Only info that I have about the original map projection is this:
"NB: the NetCDF files above use a Polar Stereographic North (EPSG:3413) projection with a horizontal resolution of 1 km x 1 km. The reference point is located at 45ºW longitude and 70ºN latitude."

is that possible to use CDO to retroject to a new netcdf file? I'm not pretty sure If I started right with remapbil or remapcon.

Any help would be appreciated

Replies (5)

RE: Project generic grid (km) to lat lon grid - Added by Helio Neto 7 months ago

So far I got some updates..

First I need to set the correct information about the projection to my netcdf file.

I ran cdo griddes smb_rec.2015-2099.BN_RACMO2.3p2-CESM2-SSP5-85.1km.YY.nc > mygrid and I got this:

#
# gridID 1
#
gridtype  = generic
gridsize  = 4039200
datatype  = float
xsize     = 1496
ysize     = 2700
xname     = x
xlongname = "x" 
xunits    = "km" 
yname     = y
ylongname = "y" 
yunits    = "km" 
xfirst    = -638956
xinc      = 1000
yfirst    = -3354596
yinc      = 1000

There is no info about the Polar Stereographic Projection (EPSG:3413)

So I add manually in the botem of mygrid file and change "generic" to "projection" in gridtype:

grid_mapping = Polar_stereographic
grid_mapping_name = polar_stereographic
proj_params = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"

And I use:

cdo -setgrid,mygrid smb_rec.2015-2099.BN_RACMO2.3p2-CESM2-SSP5-85.1km.YY.nc out.nc

And seems like everything went fine:

cdo    setgrid: Processed 343332000 values from 1 variable over 85 timesteps [5.96s 31MB]

Now my sinfo is something like this:

   File format : NetCDF
    -1 : Institut Source   T Steptype Levels Num    Points Num Dtype : Parameter ID
     1 : unknown  unknown  v instant       1   1   4039200   1  F32  : -1            
   Grid coordinates :
     1 : projection               : points=4039200 (1496x2700)
                          mapping : polar_stereographic
                      proj_params : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs
                                x : -638956 to 856044 by 1000 km
                                y : -3354596 to -655596 by 1000 km
   Vertical coordinates :
     1 : surface                  : levels=1
   Time coordinate :
                             time : 85 steps
     RefTime =  2015-01-01 00:00:00  Units = days  Calendar = standard  Bounds = true
  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss
  2015-06-30 18:00:00  2016-06-30 18:00:00  2017-06-30 18:00:00  2018-06-30 18:00:00
  2019-06-30 18:00:00  2020-06-30 18:00:00  2021-07-02 00:00:00  2022-06-30 18:00:00
  2023-06-30 18:00:00  2024-06-30 18:00:00  2025-06-30 18:00:00  2026-06-30 18:00:00
  2027-06-30 18:00:00  2028-06-30 18:00:00  2029-06-30 18:00:00  2030-06-30 18:00:00
  2031-06-30 18:00:00  2032-06-30 18:00:00  2033-06-30 18:00:00  2034-06-30 18:00:00
  2035-06-30 18:00:00  2036-06-30 18:00:00  2037-06-30 18:00:00  2038-06-30 18:00:00
  2039-06-30 18:00:00  2040-06-30 18:00:00  2041-06-30 18:00:00  2042-06-30 18:00:00
  2043-06-30 18:00:00  2044-06-30 18:00:00  2045-06-30 18:00:00  2046-06-30 18:00:00
  2047-06-30 18:00:00  2048-06-30 18:00:00  2049-06-30 18:00:00  2050-06-30 18:00:00
  2051-06-30 18:00:00  2052-06-30 18:00:00  2053-06-30 18:00:00  2054-06-30 18:00:00
  2055-06-30 18:00:00  2056-06-30 18:00:00  2057-06-30 18:00:00  2058-06-30 18:00:00
  2059-06-30 18:00:00  2060-06-30 18:00:00  2061-06-30 18:00:00  2062-06-30 18:00:00
  2063-06-30 18:00:00  2064-06-30 18:00:00  2065-06-30 18:00:00  2066-06-30 18:00:00
  2067-06-30 18:00:00  2068-06-30 18:00:00  2069-06-30 18:00:00  2070-06-30 18:00:00
  2071-06-30 18:00:00  2072-06-30 18:00:00  2073-06-30 18:00:00  2074-06-30 18:00:00
  2075-06-30 18:00:00  2076-06-30 18:00:00  2077-06-30 18:00:00  2078-06-30 18:00:00
  2079-06-30 18:00:00  2080-06-30 18:00:00  2081-06-30 18:00:00  2082-06-30 18:00:00
  2083-06-30 18:00:00  2084-06-30 18:00:00  2085-06-30 18:00:00  2086-06-30 18:00:00
  2087-06-30 18:00:00  2088-06-30 18:00:00  2089-06-30 18:00:00  2090-06-30 18:00:00
  2091-06-30 18:00:00  2092-06-30 18:00:00  2093-06-30 18:00:00  2094-06-30 18:00:00
  2095-06-30 18:00:00  2096-06-30 18:00:00  2097-06-30 18:00:00  2098-06-30 18:00:00
  2099-06-30 18:00:00
cdo    sinfo: Processed 1 variable over 85 timesteps [0.01s 15MB]

Now, when I tried to remap i got this:

cdo remapbil,r360x180 out.nc new.nc

cdo    remapbil (Abort): proj library support not compiled in!

I dont have any info about lon and lat yet.. My guess is cdo does not recgonizae "x and y" so simple to make de conversion from EPSG3413 to EPSG4326.

Any ideas?

RE: Project generic grid (km) to lat lon grid - Added by Helio Neto 7 months ago

Just some few updates, now I am capable to -setgrid correctly my original .nc file but still, I am not capable to remapbil to r360x180.

I manually set mygrid (using cdo griddes) from this:

#
# gridID 1
#
gridtype  = generic
gridsize  = 4039200
datatype  = float
xsize     = 1496
ysize     = 2700
xname     = x
xlongname = "x" 
xunits    = "km" 
yname     = y
ylongname = "y" 
yunits    = "km" 
xfirst    = -638956
xinc      = 1000
yfirst    = -3354596
yinc      = 1000

to this:

gridtype  = projection
gridsize  = 4039200
xsize     = 1496
ysize     = 2700
xname     = x
xlongname = "x" 
xunits    = "m" 
yname     = y
ylongname = "y" 
yunits    = "m" 
xfirst    = -638956
xinc      = 1000
yfirst    = -3354596
yinc      = 1000
grid_mapping = Polar_Stereographic
grid_mapping_name = polar_stereographic
standard_parallel = 60.
straight_vertical_longitude_from_pole = 315.
latitude_of_projection_origin = 90.
proj_params = "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-60 +x_0=0 +y_0=0 +k=1 +a=6378137 +b=6356257 +units=m +no_defs"

but when I try to do remapbil I got this:


cdo remapbil,r360x160 new.nc out.nc

cdo    remapbil (Abort): proj library support not compiled in!

this is the sinfo of my new.nc:

   File format : NetCDF
    -1 : Institut Source   T Steptype Levels Num    Points Num Dtype : Parameter ID
     1 : unknown  unknown  v instant       1   1   4039200   1  F32  : -1            
   Grid coordinates :
     1 : projection               : points=4039200 (1496x2700)
                          mapping : polar_stereographic
                      proj_params : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-60 +x_0=0 +y_0=0 +k=1 +a=6378137 +b=6356257 +units=m +no_defs
                                x : -638956 to 856044 by 1000 m
                                y : -3354596 to -655596 by 1000 m
   Vertical coordinates :
     1 : surface                  : levels=1
   Time coordinate :
                             time : 85 steps
     RefTime =  2015-01-01 00:00:00  Units = days  Calendar = standard  Bounds = true
  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss
  2015-06-30 18:00:00  2016-06-30 18:00:00  2017-06-30 18:00:00  2018-06-30 18:00:00
  2019-06-30 18:00:00  2020-06-30 18:00:00  2021-07-02 00:00:00  2022-06-30 18:00:00
  2023-06-30 18:00:00  2024-06-30 18:00:00  2025-06-30 18:00:00  2026-06-30 18:00:00
  2027-06-30 18:00:00  2028-06-30 18:00:00  2029-06-30 18:00:00  2030-06-30 18:00:00
  2031-06-30 18:00:00  2032-06-30 18:00:00  2033-06-30 18:00:00  2034-06-30 18:00:00
  2035-06-30 18:00:00  2036-06-30 18:00:00  2037-06-30 18:00:00  2038-06-30 18:00:00
  2039-06-30 18:00:00  2040-06-30 18:00:00  2041-06-30 18:00:00  2042-06-30 18:00:00
  2043-06-30 18:00:00  2044-06-30 18:00:00  2045-06-30 18:00:00  2046-06-30 18:00:00
  2047-06-30 18:00:00  2048-06-30 18:00:00  2049-06-30 18:00:00  2050-06-30 18:00:00
  2051-06-30 18:00:00  2052-06-30 18:00:00  2053-06-30 18:00:00  2054-06-30 18:00:00
  2055-06-30 18:00:00  2056-06-30 18:00:00  2057-06-30 18:00:00  2058-06-30 18:00:00
  2059-06-30 18:00:00  2060-06-30 18:00:00  2061-06-30 18:00:00  2062-06-30 18:00:00
  2063-06-30 18:00:00  2064-06-30 18:00:00  2065-06-30 18:00:00  2066-06-30 18:00:00
  2067-06-30 18:00:00  2068-06-30 18:00:00  2069-06-30 18:00:00  2070-06-30 18:00:00
  2071-06-30 18:00:00  2072-06-30 18:00:00  2073-06-30 18:00:00  2074-06-30 18:00:00
  2075-06-30 18:00:00  2076-06-30 18:00:00  2077-06-30 18:00:00  2078-06-30 18:00:00
  2079-06-30 18:00:00  2080-06-30 18:00:00  2081-06-30 18:00:00  2082-06-30 18:00:00
  2083-06-30 18:00:00  2084-06-30 18:00:00  2085-06-30 18:00:00  2086-06-30 18:00:00
  2087-06-30 18:00:00  2088-06-30 18:00:00  2089-06-30 18:00:00  2090-06-30 18:00:00
  2091-06-30 18:00:00  2092-06-30 18:00:00  2093-06-30 18:00:00  2094-06-30 18:00:00
  2095-06-30 18:00:00  2096-06-30 18:00:00  2097-06-30 18:00:00  2098-06-30 18:00:00
  2099-06-30 18:00:00
cdo    sinfo: Processed 1 variable over 85 timesteps [0.01s 20MB]

Any ideas?

RE: Project generic grid (km) to lat lon grid - Added by Uwe Schulzweida 7 months ago

The PROJ library is used in CDO for the coordinate transformation. Your CDO version seems to have been built without the PROJ library. When configuring/building CDO, the PROJ library is added as follows:

./configure --with-proj=<PATH_TO_PROJ_INSTALLATION> ....

RE: Project generic grid (km) to lat lon grid - Added by Helio Neto 7 months ago

This is my cdo -V:

Climate Data Operators version 2.3.0 (https://mpimet.mpg.de/cdo)
System: x86_64-apple-darwin19.6.0
CXX Compiler: clang++ -std=gnu++17 -g -O2  -pthread
CXX version : Apple clang version 12.0.0 (clang-1200.0.32.29)
C Compiler: clang -g -O2  -pthread -pthread
C version : Apple clang version 12.0.0 (clang-1200.0.32.29)
F77 Compiler: gfortran -g -O2
F77 version : GNU Fortran (Homebrew GCC 13.2.0) 13.2.0
Features: 32GB 8threads c++17 Fortran pthreads HDF5 NC4/HDF5 OPeNDAP sz sse4_2
Libraries: yac/3.0.1 NetCDF/4.9.2 HDF5/1.14.3(h1.14.2)
CDI data types: SizeType=size_t
CDI file types: srv ext ieg grb1 grb2 nc1 nc2 nc4 nc4c nc5 nczarr 
     CDI library version : 2.3.0
 cgribex library version : 2.1.1
 ecCodes library version : 2.35.0
  NetCDF library version : 4.9.2 of Mar 17 2024 23:48:59 $
    exse library version : 1.5.0
    FILE library version : 1.9.1

I installed using brew install cdo.

Can I just point the path for Proj library? Or I need to reinstall?

RE: Project generic grid (km) to lat lon grid - Added by Uwe Schulzweida 7 months ago

A new installation is required for this, configuration/build/install.

(1-5/5)

Project

General

Profile

CDO