Slight shift on remapping Polar Stereographic maps

Added by Bimochan Niraula over 3 years ago

Hi everyone,

I am trying to remap ice-concentration maps from the US Natl Ice center into OSI SAF grid. They're in GRIB format, so I simply do this:

cdo -f nc remapbil,OSI450xmpl.nc NIC.IMS_ice_v3_201433700_4km.GRIB2 Remapbil_NIC.IMS_ice_v3_201433700_4km.nc

but the resulting image has a slight offset, which can be seen most clearly near the westcoast of greenland and norway.

I also tested with other targetgrids and the offset is always there. I can also download the lat and lon data from the NIC website in their own files and when I try to plot with them, the offset is not there:

So the problem should have occured in the remapping step. But I'm not sure how exactly to solve it.

Download all files

IMS4k_directmap.png (3.14 MB) IMS4k_directmap.png
IMS4k_RemapBilOSI_fromGrib.png (1010 KB) IMS4k_RemapBilOSI_fromGrib.png

Replies (10)

RE: Slight shift on remapping Polar Stereographic maps - Added by Bimochan Niraula over 3 years ago

I also tried to creat a grid description file and remap using it. So it should be something like:

cdo -f nc remapbil,OSI450xmpl.nc -setgrid,GridDescNIC4k.txt NIC.IMS_ice_v3_201433700_4km.GRIB2 Remapbil_NIC.IMS_ice_v3_201433700_4km.nc

but I am having trouble making the grid description file. Based on the documentation of the original data* this is what I have so far:

gridtype = projection
gridsize = 37748736
xsize = 6144
ysize = 6144
xname = x
xlongname = "projection_x_coordinate"
xunits = "m"
yname = y
ylongname = "projection_y_coordinate"
yunits = "m"
xfirst = -0.5
xinc = 4000
yfirst = -0.5
yinc = 4000
grid_mapping = Polar_Stereographic_Ellipsoid
grid_mapping_name = Polar_Stereographic
longitude_of_projection_origin = 280.
latitude_of_projection_origin = 90.
semi_major_axis = 6356752.314245

but it says Unsupported projection coordinates, so clearly something is wrong.

Original documentation - https://nsidc.org/sites/nsidc.org/files/G02156-V001-UserGuide.pdf#%5B%7B%22num%22%3A31%2C%22gen%22%3A0%7D%2C%7B%22name%22%3A%22XYZ%22%7D%2C84%2C245%2C0%5D

RE: Slight shift on remapping Polar Stereographic maps - Added by Karin Meier-Fleischer over 3 years ago

Hi Bimochan,

can you upload the OSI450xmpl.nc and the NIC.IMS_ice_v3_201433700_4km.GRIB2?

-Karin

RE: Slight shift on remapping Polar Stereographic maps - Added by Bimochan Niraula over 3 years ago

Hi Karin,

I have uploaded the OSI450xmpl.nc here.
The Grib file is too large for the forum but it is available here: https://usicecenter.gov/File/DownloadProduct?products=%2Fims%2Fims_v3%2Fimsgrib2%2F4km%2F2015&fName=NIC.IMS_ice_v3_201502000_4km.GRIB2.gz

-Bimo

OSI450xmpl.nc (5.4 MB) OSI450xmpl.nc

RE: Slight shift on remapping Polar Stereographic maps - Added by Karin Meier-Fleischer over 3 years ago

Sorry, but I'm not able to get rid of the shift in the grid. Maybe Uwe can have a look at it.

RE: Slight shift on remapping Polar Stereographic maps - Added by Bimochan Niraula over 3 years ago

Nevertheless, thank you for looking into it!

I also found a maskfile with lat,lon, cellarea for the same grid : ftp://sidads.colorado.edu/DATASETS/NOAA/G02186/ancillary/masiemask_ims4km.nc (this file is ~500 MB in size). Taking the lat-lon from here and plotting the original data seems to work fine but didn't improve on the shift after remapping.

RE: Slight shift on remapping Polar Stereographic maps - Added by Uwe Schulzweida over 3 years ago

Hi Bimochan,

According to your original documentation (Section 1.8.2), the projection parameters are as follows:

+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6356257 +units=m +no_defs

These projection parameter can't be stored in GRIB format. GRIB has some restrictions to describe projections.
The above radii a=6378137 and b=6356257 are not available in GRIB. Also x_0 and y_0 cannot be stored directly,
they are stored in degrees (longitudeOfFirstGridPointInDegrees=235 and latitudeOfFirstGridPointInDegrees=-21.492).

The following output is the grid description as interpreted by CDO.

cdo griddes NIC.IMS_ice_v3_201502000_4km.GRIB2:
gridtype  = projection
gridsize  = 37748736
xsize     = 6144
ysize     = 6144
xname     = x
xunits    = "m" 
yname     = y
yunits    = "m" 
xfirst    = 0
xinc      = 4000
yfirst    = 0
yinc      = 4000
grid_mapping = Polar_Stereographic
grid_mapping_name = polar_stereographic
standard_parallel = 60.
straight_vertical_longitude_from_pole = 280.
latitude_of_projection_origin = 90.
earth_radius = 6378160.
false_easting = 12358475.909349
false_northing = 12358475.909349
longitudeOfFirstGridPointInDegrees = 235.
latitudeOfFirstGridPointInDegrees = -21.492

The earth_radius is wrong because GRIB can't store the correct radii a=6378137 and b=6356257.
false_easting and false_northing (this corresponds to x_0 and y_0) are recomputed by CDO from
the radius and longitudeOfFirstGridPointInDegrees/latitudeOfFirstGridPointInDegrees.

The correct values are calculated as follows:

proj +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +k=1 +a=6378137 +b=6356257 +units=m +no_defs << EOR
235 -21.492
EOR

result: -12286016.39    -12286016.39

false_easting and false_northing (x_0/y_0) are the absolute value of the above result!

The correct projection parameters for your GRIB file are then:

+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +x_0=12286016.39 +y_0=12286016.39 +k=1 +a=6378137 +b=6356257 +units=m +no_defs

You must put these parameters as variable proj_params at the end of the grid description. Then the geographic coordinates are calculated from them and all other parameters are ignored:

gridtype  = projection
gridsize  = 37748736
xsize     = 6144
ysize     = 6144
xname     = x
xunits    = "m" 
yname     = y
yunits    = "m" 
xfirst    = 0
xinc      = 4000
yfirst    = 0
yinc      = 4000
grid_mapping = Polar_Stereographic
grid_mapping_name = polar_stereographic
standard_parallel = 60.
straight_vertical_longitude_from_pole = 280.
latitude_of_projection_origin = 90.
earth_radius = 6378160.
false_easting = 12358475.909349
false_northing = 12358475.909349
longitudeOfFirstGridPointInDegrees = 235.
latitudeOfFirstGridPointInDegrees = -21.492
proj_params = "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +x_0=12286016.39 +y_0=12286016.39 +k=1 +a=6378137 +b=6356257 +units=m +no_defs"

If all unnecessary keys are removed it looks like this:

gridtype  = projection
gridsize  = 37748736
xsize     = 6144
ysize     = 6144
xname     = x
xunits    = "m" 
yname     = y
yunits    = "m" 
xfirst    = 0
xinc      = 4000
yfirst    = 0
yinc      = 4000
grid_mapping = Polar_Stereographic
grid_mapping_name = polar_stereographic
proj_params = "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +x_0=12286016.39 +y_0=12286016.39 +k=1 +a=6378137 +b=6356257 +units=m +no_defs"

With this grid description you can then interpolate:

cdo -f nc remapbil,OSI450xmpl.nc -setgrid,GridDescNIC4k.txt NIC.IMS_ice_v3_201433700_4km.GRIB2 Remapbil_NIC.IMS_ice_v3_201433700_4km.nc

Cheers,
Uwe

RE: Slight shift on remapping Polar Stereographic maps - Added by Bimochan Niraula over 3 years ago

Hi Uwe,

Thank you for taking the time to look into this and helping me out. I tried with the grid description you gave. It still asks for the polar_stereographic mapping parameters, so I am not sure it is using proj_params. So I copied the whole description and tried remapping, but unfortunately it still seems to have the shift.

New_Mar8First.png (1.04 MB) New_Mar8First.png

RE: Slight shift on remapping Polar Stereographic maps - Added by Uwe Schulzweida over 3 years ago

Which CDO version do you use (cdo -V)?
The grid description key proj_params is available since version 1.9.8.
Try to use the key proj4_params for older releases.

RE: Slight shift on remapping Polar Stereographic maps - Added by Bimochan Niraula over 3 years ago

Yes!! That worked.
I'm on version 1.9.6.

Thank you so much again! I have been stuck at this problem for days and finally can move on

Cheers,
Bimo

RE: Slight shift on remapping Polar Stereographic maps - Added by Bimochan Niraula over 3 years ago

Hi Uwe,

I am once again facing a similar problem with a different dataset. I am trying to convert a file from polarstereo (NSIDC) to the ease2 /lambert_azimuthal_equal_area grid (OSI SAF) grid, but the result looks very wrong. Would you mind taking a look at my steps and telling me where I made a mistake?

The user guide to the original data (https://nsidc.org/sites/nsidc.org/files/G02202-V004-UserGuide.pdf) gives this information:

PROJ4 string : +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

The spatial coordinates for the Northern polar region are the following:
Northernmost Latitude: 31.10° N
Southernmost Latitude: 89.84° N
Easternmost Longitude: 180° E
Westernmost Longitude: 180° W

Based on this, I tried to re-compute x_0 and y_0

proj +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +a=6378273 +b=6356889.449 +units=m +no_defs << EOR
> 180 31.10
> EOR
-4940045.66    4940045.66

and created this grid description file:

gridtype = projection
gridsize = 136192
xsize = 304
ysize = 448
xname = x
xlongname = "x coordinate of projection"
xunits = "m"
yname = y
ylongname = "y coordinate of projection"
yunits = "m"
xfirst = -3850000
xinc = 25000
yfirst = 5850000
yinc = -25000
grid_mapping = polar_stereographic
grid_mapping_name = polar_stereographic
proj4_params= "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=4940045.66 +y_0=4940045.66 +a=6378273 +b=6356889.449 +units=m +no_defs"

With this I tried to interpolate:
cdo -f nc remapnn,OSI450xmpl.nc -setgrid,CDRv4GrdDescNH.txt -selvar,cdr_seaice_conc NSIDC_20200101_f17_v04r00.nc remapNN_NSIDC_20200101_f17_v04r00.nc

but the result is scrambled so I'm obviously doing something wrong.

I have attached the sample files OSI450xmpl.nc and NSIDC_20200101_f17_v04r00.nc here.

Download all files

OSI450xmpl.nc (9.4 MB) OSI450xmpl.nc
NSIDC_20200101_f17_v04r00.nc (227 KB) NSIDC_20200101_f17_v04r00.nc

(1-10/10)

Project

General

Profile

CDO